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ABSTRACT 

As two black holes bound to each other in a close binary approach merger their inspi- 
ral time eventually becomes shorter than the characteristic inflow time of surrounding 
orbiting matter. Using an innovative technique in which we represent the changing 
spacetime in the region occupied by the orbiting matter with a 2.5PN approximation 
and the binary orbital evolution with 3.5PN, we have simulated the MHD evolution 
of a circumbinary disk surrounding an equal-mass non-spinning binary. Prior to the 
beginning of the inspiral, the structure of the circumbinary disk is predicted well by 
extrapolation from Newtonian results. The binary opens a low-density gap whose ra- 
dius is roughly two binary separations, and matter piles up at the outer edge of this 
gap as inflow is retarded by torques exerted by the binary; nonetheless, the accretion 
rate is diminished relative to its value at larger radius by only about a factor of 2. 
During inspiral, the inner edge of the disk at first moves inward in coordination with 
the shrinking binary, but as the orbital evolution accelerates, the rate at which the 
inner edge moves toward smaller radii falls behind the rate of binary compression. In 
this stage, the rate of angular momentum transfer from the binary to the disk slows 
substantially, but the net accretion rate decreases by only 10-20%. When the binary 
separation is tens of gravitational radii, the rest-mass efficiency of disk radiation is a 
few percent, suggesting that supermassive binary black holes in galactic nuclei could be 
very luminous at this stage of their evolution. If the luminosity were optically thin, it 
would be modulated at a frequency that is a beat between the orbital frequency of the 
disk's surface density maximum and the binary orbital frequency. However, a disk with 
sufficient surface density to be luminous should also be optically thick; as a result, the 
periodic modulation may be suppressed. 



Subject headings: Black hole physics - magnetohydrodynamics - accretion, accretion 
disks - Galaxies: nuclei 
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Introduction 



There is now excellent evidence tha t every galaxy with a bulge contains a supermassive black 
hole at its center (jGiiltekin et al.ll2009l ). In addition, the prevailing theory of galaxy formation 
posits that today's massive galaxi es were assembled from smaller pieces, as dark-matter haloes of 
progressively greater size merged (IDavis et al.l Il985l : iBardeen et al.l Il986l ) . If massive black holes 
were already present in those progenitors, they would bring their black holes with them into the 
new combined galaxy, creating an opportunity for the black holes to merge. Such an event would be 
very exciting to detect for many reasons: It would reveal the presence of supermassive black holes 
early in the life of galaxies. It would shed important li ght on the growth of the strong correlations 
between nuclear black hole mass and galaxy structure (|Gurtekin et al.ll2009l ). Most of all, it would 
provide a concrete example of one of general relativity's most spectacular predictions and possibly 
also allow a test of the validity of general relativity in a truly strong-field regime. 

An extremely large amount of energy is very rapidly released in a binary black hole (BBH) 
merger event, almost all of it through gravitational radiation [several percent of the black hole 
masses in a timescale of ~ (Mbbr/Mq) x 493 fjs]. Gravitational radiation may be strong enough 



to eject the final remnant from its host galaxy, w ith recoil ve 



km/s predicted bv numerical relativitv simulations (Baker et al. 


2007; 


CamDanelli et al. 


2007a 


b; 


Gonzalez et al. 


200' 


7. 


Herrmann et al. 




2007 


; ] 


<!oDDitz et al. 


2007 


Baker et al. 


2008; 


Healv et al. 


2009 


■ 


Lousto et al. 


2010; 


Lousto & Zlochowei 


2011: 


Lousto et al. 


201^ 


). Unfortunately, a gravi- 



ociti e s or "kicks" up to 1 3 



tational wave observatory with adequate sensitivity in the appropriate frequency r ange is still well 
in the future, w hether it operates by direct detection or through pulsar timing (|Jennrichl |2009| ; 
Wen et al.ll201ll ). On the other hand, even if only a small part of the energy is deposited in nearby 



gas, the associated photon signals might be much more readily seen with instruments operating 
today. Because the energy given to the gas comes from work done by gravitational forces, one would 
expect, on the basis of the Equivalence Principle, that the total energy added to the gas would be 
proportional to its mass. If most of this added energy is dissipated into heat (local irregularities are 
likely to drive sh ocks), the to tal energy radiated in photons would then be similarly proportional 
to the gas mass (JKrolikJ |2010j) . The question is, therefore: "How much mass would one expect in 
the neighborhood of a black hole merger?" 

Even if a BBH were supplied with mass at a rate characteristic of high luminosity quasars (~ 
IOMq yr _1 ), several effects may severely reduce how much gas remains close to the binary. Torques 
exerted by the binary on the inflowing gas may hold back the inflow, preventing much of it from ap- 



proa ching closer than a few times the binary separation a (jPringle 



1991 



MacFadven Milosavljevic 



2008). As the binary compresses, whether by interactions with passing stars and external gas or 
by gravitational radiation, the gas follows, but is held off at a distance of at least ~ 2a. Toward 
the end of the binary's evolution, gravitational radiation losses grow rapidly and dominate the 
orbital shrinkage. Ultimately, the orbit shrinks on a timescale shorter than the characteristic ac- 
cretion inflow time and the BBH is expected to decouple from the disk. After such decoupling, 
there would not be enough time for much disk mass to catch up with the black holes before they 
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merge (|Milosavlievic k Phinnevll2005l ). Thus, for a given external supply rate, the amount of gas 
available to be heated in a merger is determined by a competition between the internal stresses 
that drive inflow and a pair of dynamical mechanisms that tend to keep gas at "arms-length" from 
the merging black holes. 

Until recently, efforts to quantify these effects have relied alm ost entirely on the phenomeno- 
logical Shakura-Sunyaev a-disk model to describe internal stresses (IM ilosavli evic k Phinnevll2005l ; 
MacFadven k Milosavljevidbood : iLiu k Shapirohoid : banaka k Menoulhoioh . in which the verticallv- 
integrated and time- and azimuthally-averaged internal stress is supposed to be a factor a times the 
similarly integrated and averaged pressure; the only exceptions were studies focusing on binaries 
with large mass ratios between primary and seco ndary, a limit primar i ly relevant to planet forma- 



tion and to extreme-mass ratio inspiral sources ( Winter s et al 



Nelson k Papaloizou 
Farris et all (|201lf 



2003 



and 



Kocsis et al. 2011: Yunesetal 
II 11 1 " 1 



2011 



20031 : iPapaloizou k Nelson 



2003 



Moreover, with the exception of 
Bode et al.l (|2012l ). which assumed these stresses were negligible, all these 



calculations also assumed Newtonian dynamics. However, there is strong reason to think that 
the actual mechanism of these stresses is magnetohydrodynam ic (MHD) turbulence, stirred by the 
magnetorotational instability (MRI) (jBalbus k Hawlevl Il998l ). In contrast to ordinary, isotropic 
turbulence, orbital shear makes this turbulence highly anisotropic, so that there is a non-zero cor- 
relation between the radial and azimuthal components of the magnetic field; this correlation creates 
the stress. MHD calculations are therefore required, at the very least, to establish the appropriate 
scale of the stresses and the approximate magnitude of a. In addition, although the a-model may 
give a reasonable description of time-averaged behavior well inside the body of an a ccretion flow. 



it is particularly ill-suit ed for predicting dynamical behavio r on shorter timescales (jHirose et al 



20091 ) and at disk edges (jKrolik et alJl2005l ; iNoble et al.ll2010l ). Because the key issues in how much 
gas reaches a merging BBH depend, of course, on the time-dependent behavior of gas near and 
within the inner edge of a disk, explicit calculation of the MHD turbulence is also required for an 
accurate treatment of the time- and spatial-dependence of the internal stress. 

In this paper, we present the first simulation of a circumbinary accretion disk around a binary 
black hole system during the epoch in which the binary's inspiral time grows shorter than the inflow 
time through the disk. Generically, this period occurs not long before the binary's final merg e r. Ou r 



physics treatment includes fully relativistic MHD. This study differs from that of lShi et al.l (|2011 



who presented similar MHD simulations, but concentrated on the Newtonian regime, when the 
black holes are very widely separated. Moreover, their Newtonian treatment did not allow for the 
black holes to inspiral, a reasonable assumption when the semi-major axis is hundreds of thousands 
of gravitational radii, but a terrible assumption in the late inspiral. We here focus on BBHs with 
separations of ~ (10-20) r g , where r g = GM/c 2 and M is the total mass of the binary (we adopt 
geometric units with G = c = 1 for the remainder of this paper). The spacetime associated with 
the BBH orbital dyn amics is d e scribe d through a vacuum post-Newtonian (PN) approximation (see 
the review paper of iBlanchetl (120021 ) and references therein), where we neglect the back reactio n 
of the disk on the BBH dynamics. Our work also contrasts with that of iGiacomazzo et al.l (J2012J), 
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who employed full numerical relativity to compute the spacetime in which an initially uniform gas 
distribution with an externally-imposed uniform magnetic field evolved during the last 3 orbits 
before merger. 

As we will describe below, the PN approximation is adequate to describe the spacetime evo- 
lution for our needs. The PN scheme is a method to describe approximately the dynamics of 
physical systems in which motions are slow compared to the speed of light and gravitational fields 
are weak. That is, one solves the Einstein field equations perturbatively, expanding in (v/c) 2 <C 1 
and r g /r = GM/(rc 2 ) <C 1 Here v, M and r are the characteristic velocity, mass and size or 
separation of the system. This approximation has been remarkably effective in describing the 



perihelion precession of Mercury lEinsteinl (|1915l ). and the gravita tional- wave loss from bina ry sys- 



tems, such as the Hulse-Taylor pulsar, PSR B1913+16 (see e.g. Iweisberg fc Tavlorl (feoosl ): Iwill 



(|201ll )) . PN theory also plays a key role in the construction of the gravitational-waveform tem- 



plates (jSathvaprakash &; Schutal2009l ) for inspiraling compact objects currently used in the search 



for gravitational waves by laser-interferometric observatories. PN theory has also been recently 



mergers ( 


Tichv et al. 


2003; 


Bonning et al. 


2003; 


Yunes et al. 


2006; 


Yuncs Sz Tichy 


2006 




Yunes 


2007 




Kellv et al. 


200^ 




CamDanelli et al. 


2009 




Johnson-McDaniel et al. 


2009; 


Kellv et al. 


2010; 


Mundim et al. 


2011). In all cases, the PN approximation is developed to sufficiently high pertur- 



bative order that the error contained in the approximation is much smaller than that associated 
with either the data in hand (in the case of binary pulsars) or the data expected (in the case of 
direct gravitational wave detection). 

Using this PN-approximated description of the spacetime, we first evolve the BBH at a fixed 
initial separation ao = 20M to allow the accretion disk to relax to a quasi-steady state. To study the 
effect of orbital shrinkage on the accretion disk, we then follow the binary inspiral until it reaches a 
separation a = 8M, beyond which the PN approximation ceases to be sufficiently accurate for our 
purposes. In a separate simulation, we kept the binary's separation fixed at 20M and continued to 
evolve until ~ 76000M to study the secular dynamics of the quasi-steady state, and, by contrasting 
with the first simulation, highlight the special effects induced by the inspiral. 

Our findings can be summarized as follows. The mass at r ~ 2.5a builds steadily throughout 
the quasi-steady state, but much of it eventually concentrates in a distinct "lump" . At smaller radii, 
a gap is cleared as torques and forces exerted by the binary either sweep matter inward or fling it 
outward. Much of the small amount of mass in this gap is found in a pair of streams emanating 
from the inner edge of the disk and curving inward toward each black hole. These streams carry 
nearly half of the mass accreting through the bulk of the disk to the inner boundary of the problem 
volume at r = 0.75ao- 

As the binary starts to shrink, the inner edge of the disk at first moves inward following the 
orbital evolution of the binary, but eventually cannot keep up, as the orbital shrinkage grows faster. 
Nevertheless, a significant amount of mass still follows the binary's inspiral within the gap region. 
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We find that the final accretion rate in the inspiral stage is about 70% of the corresponding rate 
in the steady-state stage. 

The luminosity of the disk is proportional to its surface density. If the accretion rate fed to the 
disk were comparable to that of ordinary AGN, the surface density and therefore the luminosity, 
of such a circumbinary disk could approach AGN level. Most strikingly, the luminosity should be 
modulated periodically at a frequency determined by the binary orbital frequency and the binary 
mass ratio — 1.46 times the binary orbital frequency in the case of equal masses. However, the 
amplitude of modulation may be reduced by the large optical depth of the disk if the surface 
density is large enough to generate a sizable luminosity. 

This paper is organized as follows. In Section [2] we describe the construction of the dynamical 
BBH spacetime. In Section [3l we report the details of the MHD simulations of the disk. In 
Section [3] and Section [5j we present the results from these simulations and interpret them in the 
light of previous work and from the point of view of potential observational signatures. Finally, in 
Section [6l we summarize our principal conclusions. 



2. Binary Black-Hole Spacetime 

While solutions of th e Einstein equations for single black holes were discovered as early as 
1916 (j Schwar zschildl 1 1 9 1 6l ) . no exact closed-form solution to the two-body problem exists and one 
generally needs to solve the Einstein equation numerically. With the breakthroughs in numerical 



relativity (jPretoriusI 120051 : ICampanelli et al.l 120061 : iBaker et al 



2006 



Scheel et al 



20061 ). it is now 



possible to perform stable and accurate full numerical simulations of BBHs in vacuum for a wide 
variety of mass ratios and spins parameters. However, because the Einstein equations can be 
thought of as modified wave equations, with wave speeds of c, the Courant condition^ greatly 
limits the timestep size, making full numerical simulations impractical when the characteristic 
MHD speeds are significantly smaller than c. On the other hand, if an approximate, but accurate, 
spacetime is given, the Courant condition is set by the MHD speeds, allowing for a much larger 
timestep. Fortunately, analytic perturbative techniques have been successfully developed to tackle 
the spacetime problem both in the regime where the black holes are not too close, as well as in 
the close limit regime where the spacetime can be treated as a perturbed single black hole. In this 
paper, we use the PN approach to model the spacetime of an inspiraling binary system prior to 
merger, neglecting the effect of the disk on the evolution of the black holes, from orbital separations 
of 20M down to ~ 8M, roughly where the standard PN approximation becomes inaccurate for our 
purposes. 



Using this PN-approximated solution, we then solve for the relativistic MHD evolution of the 



The Courant condition is a stability condition relating the timestep dt to the spatial resolution h. The timestep 
is limited by dt < h/v, where v is the fastest propagation speed in the system of interest. 
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circumbinary accretion disk. We stop the PN evolution at r = 8M, but one can in principle 
continue the simulations beyond this regime. To do so, one could use a snapshot of the PN metric 
and MHD data at that radius as initial conditions to then carry out a fully non-linear GRMHD 
evolution, using numerical relativity techniques to solve the coupled GRMHD Einstein system of 
equations (which will be the subject of an upcoming paper). 

We perform two simulations: (i) RunSS keeps the semi-major axis of the binary artificially 
fixed at 20M; (ii) Runln starts from a snapshot of RunSS at t = 40000 M (or after ~ 70 orbits) 
and then lets the black holes inspiral at the PN-theory prescribed rate down to a separation of 
~ 8M. RunSS is used to study the secular evolution of the accretion disk at fixed binary separation, 
while Runln is used to investigate how the diminishing separation alters this secular evolution. We 
describe our PN approach to model the spacetime metric below. 



2.1. The post-Newtonian Approximation 



The PN approximation is based on a perturbative expansion of all fields, assuming slow mo- 
tion v/c -C 1 and weak fields GM/(rc 2 ) <C JE These assumptions allow us to search for solutions 
that can be expressed as a divergent asymptotic series about a flat Minkowski background space- 
time. These perturbations obey differential equations determined by the PN-expanded Einstein 
field equations. One then solves such equations perturbatively and iteratively to construct an 
approximate solution. 

The two body problem in the slow- motion/ weak-field limit is better understood by classifying 
the spacetime into di fferent regions , wher e different as s umptions hold and d ifferent approximation s 
can be used ( see e.g.lThornel (ll98oh: lAlvJ (bood . booah : lYunes et all 120061) : lYunes Tichvl 120061): 
Yunesl (|2007l ) ; iJohnson-McDaniel et al.l (|2009l ) for a review) . Here we concentrate on the near zone , 
which is the region sufficiently far from the horizons that the weak-field approximation of PN is 
valid, but less than a reduced gravitational wave wavelength A away from the center of mass of the 
system, so that retardation effects can be treated perturbatively. We note that in the far zone, i.e. 
the radiation zone where retardation effects can no longer be treated perturbatively, a multipolar 
post-Minkowskian expansion can be used rather than a PN one. Very close to each black hole 
(i.e. in the inner zone), perturbed Schwarzschild solutions are used (which can be extended to 
include spin by using perturbed Kerr solutions). In this paper, the binary black hole metric will 
be approximated with only the near zone solution. 

The PN approximation, of course, has its limits: binary systems eventually become so closely 
separated that a slow- motion/weak- field description is inappropriate. For example, consider the 
simple case of a test particle spiraling into a non-spinning (Schwarzschild) black hole in a quasi- 
circular orbit. Eventually, the particle will reach the innermost stable circular orbit, at which point 



2 Note that we have explicitly re-introduced c and G in this section in order to discuss the PN approximation. 
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its orbital velocity u/c ~ 0.41. Clearly, such a velocity is not much less than unity, and thus, 
the PN approximation need not be an accurate description of the relativistic orbital dynamics. 
Similarly, when comparable-mass BBHs inspiral, they eventually reach a separation at which the PN 
approximation is a bad predictor of the dynamics, since the small- velocity /weak- field assumptions 
are violated. 

The determination of the formal region of validity of the PN approximation is crucial, but it can 
only be assessed when one possesses a more accurate, perhaps numerical, description of the orbital 
dynamics. This is indeed the case when considering extreme mass-ratio inspirals (EMRIs), consist- 
ing of a stellar-mass compact object spiraling into a supermassive black hole. When considering 
EMRIs, one can model the spacetime through black hole perturbation theory, i.e., by decomposing 
the full metric as that of the supermassive black hole plus a perturbation in duced by the s tellar- 
mass compact object, without a ssuming slow motions (s e e e.g. Section 4 of Hughes ( 20 091) fo r a 



recent review, Mino et al.1 (|l997h ; ISasaki k Tagoshil ((2003) ; Baracki ((20081 ) ; IPoisson et al.l f)201lh for 



related topics). To leading-order, the orbital dynamics are then described by geodesies of the small 
object in the spacetime of the supermassive black hole. The orbital motions are slowly perturbed 
by the radiation reaction due to the emission of gravitational waves. 

The comparison of black hole perturbation theory and PN theory predictions has allowed for 
the construction of different measures to estimate the PN region of validity. When considering the 
reduction in signal-to-noise ratio induced by filtering an "e xact" black-hole perturbation theory 
gravitational wave with a order filter, IPoisson! (jl995l ) found that PN theory is sufficiently 

accurate provided v <^ 0.2. When considering the 5.5 and 4PN order predictions for the loss of the 
binary's binding energy for non-spinning and spinning bac kground black h o les re specti vely, relative 



to an "exact" black hole perturbation theory prediction, mines &; Bertil (|2008l ) and IZhang et al 



((201 ll ) found that the former is accurate provided v <^ 0.29, which corresponds to an orbital 
separation of a ;> 11M. The difference between these estimates is due to the different measures 
used and the different order of the PN approximation employee^. 

The region of validity of the PN approximation for comparable-mass binaries has not been as 
well studied. This is because black hole perturbation theory is not applicable here, and one must rely 
on full numerical relativistic simulations. Currently, state-of-the-art simulations can only model the 
last few tens of orbits prior to merger, while the determination of the formal region of validity would 
require knowledge of at least the last thousand orbits. Nonetheless, there exist analytical arguments 
suggesting that the region of validity in the comparable-mas s case is larger than in the EMRI 



case (i.e., the PN e xpansion is valid for even larger velocities) ((Simone et al.l 119971 ; iBlanchet 



Mora & Willi 12004 ). Moreover, the NINJA (Numerical INJection Analysis)-2 project (jAjith et al 



2003: 



3 This was extended to 5.5PN order in the erratum and addendum of iPoissonl ( 199a ). 

4 It is not surprising that beyond 3PN order the region of validity of the PN approximation shr inks. This is a 
prope rty o f divergent asympto tic series, whose behavior in the context of PN theory was analyzed bv lYunes fc Berti 
(2008) and lZhang et al.1 fcoilh . 
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20121 ). a collaboration between numerical relativists and gravitational wave data analysts, has 



established that certain 3PN order gravitational waveforms are sufficiently accurate for use as 
templates provided v <; 0.33 (a ^ 8M). 



2.2. Near-Zone PN Evolution 



The PN order of a given term is determined by the exponent of the perturbation parameter 
contained in that term. In the near zone, the PN expansion of a BBH spacetime metric is a series 
expansion in the orbital velocity v/c <C 1 and the field strength GM^/^r^c 2 ), where Ma and r& are 
the masses of the A th particle and the distances from the A th particle to a field point respectively. 
Here, we may consider 1/c as the PN parameter which goes to zero in the Newtonian limit c — > oo. 
Notice that by the virial theorem v 2 /c 2 = 0[GM/(ac 2 )]. A term proportional to (l/c) n beyond 
the Newtonian (leading-order) expression is said to be of (n/2) th PN order. 

The near zone metric will be described here by a resummed PN expression. One begins 
with the 2.5PN expansion of the metric for non -spinning point partic les in a quasi-circular orbit 
in harmonic coordinates, given for example in iBlanchet et al.1 (119981 ). Such a metric, however, 
descr ibes black holes a s poin t -particles, wh ic h is why one then applies a "bac kground resummation" , 



as m 



Yunes & Tichvi (|2006l ); lYunesI (|2007l ); IJohnson-McDaniel et al.l (|2009l ). This resummation is 



intended to improve the strong- field behavior of the metric close to each point-particle, i.e., it 
recovers the horizon of each individual black hole. The metric can then be formally written as 



g^u(t, x) = g^ u [x; yA(t), v A (t)] , 



(1) 



where x is a spatial vector from the binary's center of mass to a field point, while yA(t) an d VA(t) 
are the particle's spatial location and 3-velocity with A = (1, 2). This metric depends on the mass 
of each individual black hole, but also on the binary orbital evolution {^(i), VA(t)} that must be 
prescribed separately. We will use Greek letters (e.g., //, v, A, k) to represent spacetime indices 
[0, 1, 2, 3], and Roman letters (e.g., i, j, k, I) to represent spatial indices [1, 2, 3]. 

The orbital evolution can also be prescribed within the PN approximation. We simplify 
the analysis by considering only quasi-circular orbits. This simplification is justified because 
gravit ational wave emission tends to circular ize binaries very efficiently, as demonstrated i n the 
weak (IPetersI 11964 ) and strong field regimes (jSperhake et al.l 120081 : iHinder et al.l 120081 . l2010l ) . We 
will here use a 3.5PN expans ion for the orbital phase evolution 4>(t), as given for example by Equa- 
tion (234) of iBlanchetl (|2002j), which depends on the quantity = v(t c — t)/(5M), where t is time, 
t c is the time of coalescence, M is the total mass and v = M1M2/M 2 is the symmetric mass ratio. 
The PN orbital frequency is calculated from fibin = d(j>/dt in this paper. 

Although waveforms can be fully characterized by harmonics of the orbital phase, the latter 
depend on the orbital trajectories. One may model these in harmonic coordinates via 



V\{t) = -j^ait) [cos0(t), sin 0] , 



(2) 
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Mi 

VW) = --^o(t)[cos0(t), sin0(i), 0] 



(3) 



where a(t) = |yi(i) — 2/2 ) | is the orbit al separation as a function of time. This separation can be 
calculated via the balance law (see e.g., (|Blanchetll2002l )). which states that the local rate of change 
of the binary's orbital binding energy is exactly balanced by the gravitational wave luminosity 
carried out to future null infinity and into black hole horizons, namely 



dt 



-C. 



The rate of change of the orbital separation is then given by 



da 
~dt 



dEp 
da 



rb 



c. 



Assuming the initial condition a(t = 0) = ao, we then find 



t = t c 



\ 1 f dE o 
da 



rb 



V da' 



C 



-1 



(4) 



(5) 



(6) 



where the integrand is expanded in a Taylor PN series and the coalescence time t c is defined by 



da ( dE ° vh 
V da 



c- 



(7) 



In this paper, we only use the part of C that is carried to future null infinity and, although 
Equations (|6|7p are typically Taylor expanded to evaluate t(a), here they are inverted through 
Newton- Raphson minimization to yield a(t). 

Given the above analysis, we can now calculate the time of coalescence and the number of 
orbits in each of the simulations carried out. As already discussed, RunSS is artificially kept at 
a fixed semi-major axis, so a(t) = ao for all times, and thus, formally t c = 00. On the other 
hand, Runln keeps a(t) fixed to a(t) = ao = 20M for t < i s hrink = 40000M, after which it is 
allowed to evolve according to the PN equations of motion. The time of coalescence for this run 
can be computed by invertin g Equation (HQ) to obtain t c ~ 14000M, although to leading order it is 
approximately described by (|Peterslll964h 



tc. 



5 /oqV 
256u \MJ 



M 



5 /a \ 4 + < 
256 \MJ q 



-M . 



(8) 



where q = M2/M1 is the binary mass ratio. A BBH clearly takes longer to merge for systems 
that start at larger initial separations and that possess extreme mass ratios. Obviously, t c is always 
defined as the length of time to coalesce, when the binary is allowed to inspiral. The total simulation 
time of Runln is then t c + ^shrink ~ 54000M. 



We have plotted a few diagnostics to get a sense of the evolution of the binary system in 
Runln . Figure Q] shows the orbital evolution of the binary in the x-y plane, after it is allowed 
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Fig. 1. — The orbital motion of one of the black holes in the binary from Runln . Its trajectory 
starts from an initial separation of a(0) = 20M and stops at a(tf) ~ 8M. Its black hole companion 
is located at a parity-symmetric point across the origin (track not drawn in the figure for the sake 
of clarity). 



to inspiral. Figure [2] plots the orbital separation as a function of time. Observe that initially the 
semi-major axis is artificially kept fixed, while after t > i s hrink it is allowed to decrease due to 
gravitational radiation reaction. Figure [3] plots the number of orbits ^Vorbits traced by the binary 
system as a function of time, which is given by 



N. 



orbits 



— ^bin,0 t (if t < ishrink) 

TT" [^bin,0 ^shrink + <j>(t 
Z7T 



(9) 



^shrink)] (if t ^ ^shrink) 



We have here defined ^bin,o = ^bin(^ = oo) to be the (constant) PN orbital frequency at a fixed 
semi-major axis, and we have set 0(0) = for the PN phase evolution. The number of orbits is 
obviously a piece- wise function since when t < t s hrinki ^orbits increases linearly, as the binary is 
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Fig. 2. — The evolution of the orbital separation with respect to time, a{t). We turn on the 
gravitational radiation reaction at t = 40000M. 

artificially kept at fixed do, while when t > t s hrink, -^orbits can be approximated to leading order by 

\2 



N, 



orbits 



1 /O0\5/2 _ J_ /Oo\5/2 (1+qY 



64vri/ \MJ 



64vr \MJ 



(10) 



Therefore, for Runln there are approximately 100 total orbits, while for RunSS there are approxi- 
mately 127 orbits. 



3. Simulation Details 

Like accretion disks around single black holes, circumbinary accretion flows are well described 
by the ideal MHD equations of motion (EOM) in the curved spacetime of only the black hole 
or holes. We therefore neglect the matter's contribution to spacetime curvature and the accu- 
mulation of mass and momentum by the black holes from gas accr e tion. Many codes have bee n 



DeVi 



hers fc Hawlev! (12 003): 



written to simulate the single black hole case fe.g..lKoide et al. (1999); 

Garnmie^~al](l2003h; [Noble et al.1 |20od ); lAnninos et al.1 ( 20051 ) ; iKomissarov j 20051); Tchekhovskoy et al 
(I2OO7I): iNoble et al.l (hood)), while only t he equations of electrodynamics (jPalenzuela et al.l 120 09). 



force-free MHD ( jPalenzuela et al.l boiObl lah and nonmagnetized hydrodynamics (e.g. 



Bode et al 



-12- 




Fig. 3. — The number of orbits as a function of time. (Black) RunSS . (Grey) Runln , in which 
the first part from t = to 40000Af shows the circular orbit without the radiation reaction. 



(12O10MFarris et al.1 (j20ld . 



2011 



Bode et al.l (|2012l )) have been solved in the relativistic circumbi- 



nary setting. Unfortunately, these latter simulations employ methods, like block-structured adap- 
tive mesh refinement (AMR) in Cartesian coordinates, that typically lead to poor conservation 
of fluid angular momentum and excessive dissipation at refinement boundaries. These two effects 
alter the disk's angular momentum transport mechanism and thermodynamics in a nontrivial way. 
Furthermore, they require the solution of the Einstein equations, which — in turn — imposes a sig- 
nificant computational burden. In order to avoid these problems, we take an alterna te route and 



solve the MHD EOM using a code designed for single black hole systems: Harm3d (jNoble et al 



20091 ). Fortunately, Harm3d was written to be almost independent of coordinate system or choice of 
spacetime, so modifying it to handle non-axisymmetric, time-depende nt spacetimes was st raightfor- 
ward. In fact, the only differences between the algorithm described in lNoble et al.l (|2009l ) and here 
are that the metric (and its affine connection or gravitational source terms) needs to be updated 
every sub-step of the second-order Runge-Kutta time-integration procedurqj. Below, we describe 
the equations solved, initial data setup and other details of the disk evolution. 



5 Note that many other technical changes were made that do not affect the algorithm, but do affect the runtime 
efficiency and design of the code. 
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3.1. MHD Evolution 



Since we assume that the gas does not self-gravitate and alter the spacetime dynamics, we 
need only solve the GRMHD equations on a specified background spacetime, g^ u {x x ), where {x x } 
represents a set of general spacetime coordinates. The EOM originate from the local conservation 
of baryon number density, the local conservation of energy, and the induction equations from 
Maxwell's equations (please see iNoble et al.l (|2009l ) for more details). They take the form of a set 
of conservation laws: 

<9 4 U (P) = -d i F i (P) + S (P) (11) 



where U is a vector of "conserved" variables, F* 
Explicitly, these are 



are the fluxes, and S is a vector of source terms. 



U(P) 
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F l (P) 
S(P) = 



pu\ T't + pu' 



pu 



T 



b l u k 



b k W 



-9 



0, T\T x tK 



Ft, T K \F j K — Fj, 



(12) 
(13) 
(14) 



where g is the determinant of the metric, T X ^ K is the metric's affine connection, B l 



'4tt is our 

magnetic field (proportional to the field measured by observers traveling normal to the spacelike 
hypersurface) , *F^ U is the Maxwell tensor, u M is the fluid's 4-velocity, b» = X {8» v + vPu v ) B v 
is the magnetic 4-vector or the magnetic field projected into the fluid's co-moving frame, and 
ut / \J ~ 9 tt is the fluid's Lorentz function. The MHD stress-energy tensor, T^, is defined as 



W 



- [IV 



(ph + 2p m ) u^Uy + (p + 

Pm) 9[iv b^b u 



(15) 



where p m = b^b^/2 is the magnetic pressure, p is the gas pressure, p is the rest-mass density, 
h = 1 + e + p/p is the specific enthalpy, and e is the specific internal energy, We evolve the quantity 
(pu* + T_l} instead of TK m order to reduce the magnitude of the internal energy's numerical 



error ( Gammie et al.l 120031 ). Note that the terms proportional to r\ K and r A ^ K in the source no 



longer vanish as the metric is now dependent on time and azimuthal coordinate, 0. Also, note 
that we add a negative source term (— J 7 ^) to the local energy conservation equation to model 
energy/momentum loss from radiative cooling; please see Section [3.41 for more details. 

The MHD evolution is facilitated by calculating and using so-called primitive variables: the 
rest-mass density (p), the internal energy density (u = pe), the velocities relative to the observer 
moving normal to the spacelike hypersurface, u l = u l — u t g u /g tt . The magnetic field B l is con- 
sidered both a primitive and a conserved variable. We employ piecewise parabolic reconstruc- 
tion of the primitive variables for calculating the local Lax-Friedrichs flux at each cell interface 
(jGammie et al.ll2003l ). We use a 3 -dimensional version of the FluxCT to impose the solenoidal con- 
straint, diyJ—gB l = (|T6thll2000l ). The EMFs (electromotive forces) are calculated midway along 
each cell edge using piecewise parabolic interpolation of the fluxes from the induction equation. A 
second-order accurate Runge-Kutta method is used to integrate the EOM using the method of lines 
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once the numerical fluxes are found. The prim itive variables are found from the conserved variables 



using the "2D" scheme of iNoble et al.l (|2006). A conservation equation for the entropy density is 
evolved and used to replace the total energy equation of the 2D method whenever the plasma 
becomes too magnetically dominated, or — specifically — when pe < 0.02p m ; this procedure h e lps u s 



avoid numerical instabilities and negative pressures from developing. Please see lNoble et al.l (|2009l ) 
for more details. 

The MHD evolution is performed in the same way as in single black hole cases except that 
the metric is evaluated^] at the present sub-step's time before the MHD fields are updated. The 
metric is required in many facets of the update procedure. For example, it is used in calculating 
the 4- velocity from the primitive velocities, source terms and geometric factors in the EOM, and for 
deriving the primitive variables from the conserved variables. The affine connection is calculated 
via finite differencing the PN metric to evaluate 

^uk = \g^ a (d u g.a + d K g va - d a g UK ) . (16) 

The spatial finite differences use fourth-order centered stencils away from the physical boundaries, 
and backward/forward stencils adjacent to the physical boundaries. Since the metric is evaluated 
and stored at the cell centers and faces, but the connection is only evaluated at the centers, fourth- 
order stencils require only three cells' worth of data to compute. The time derivatives are second- 
order accurate, but use a time spacing 10 -3 times that used in the MHD integration. This means 
that additional evaluations of the metric are made at advanced and retarded times at the cell 
centers to calculate the time derivatives for each connection evaluation. We have verified that 
the truncation error from the time derivatives is smaller than that from the spatial derivatives. 
Also, the connection's spatial finite differencing is one order more accurate than that of the MHD 
procedure, implying it is not the primary source of error in the calculation. Please see Appendix IB1 
for a discussion on our resolution tests. 



3.2. Initial Conditions 



In this project, we avoid evolving the gas in the neighborhood of the black holes, choosing 
instead to focus on establishing reasonable prior conditions for the gas that ultimately feeds the 
BBH. We therefore excise a spherical domain, which includes the binary, from our calculation. 

It is common to begin with initial conditions devoid of large transient artifacts. This is often 
done by starting from a torus of material in equilibri um (via pressure and rotational support ) 
about the central gravitating source (e.g., a black hole) (|Chakrabartilll985l ; iDe Villiers et al.l l2003). 
Unfortunately, such tori will not be near equilibrium in our spacetime as it is (t, ^-dependent. 



6 We remind the reader that the metric is known in closed form, requiring only direct evaluation except for the 
Newton-Raphson iteration to find the current time's binary separation. 
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Plus, the equations describing their structure assume that the metric has the same form (i.e. share 
the same zero- valued elements) as the Kerr metric in Boyer-Lindquist coordinates. We resolve these 
issues in the following way. First, since we hold the binary at fixed separation for several orbits, the 
spacetime initially has a helical Killing symmetry, with Killing vector KL a = f2bin {d<j>) a + (dt) a '■ In 
other words, the spacetime is invariant in a frame rotating with the binary, while the separation is 
held constant. Since the torus will lie a few ao away from the binary, its dynamical response time — 
comparable to its orbital period — will be longer than the binary period, implying that a torus near 
equilibrium in this helically-symmetric spacetime will also be near equilibrium in its time average. 
Due to its helical symmetry, its time average is also its azimuthal average. We therefore start with 
a torus in equilibrium in a (t, (^-independent spacetime, g^ u , found by averaging over (f>: 



(17) 



We have verified that the same components that are zero-valued in Boyer-Lindquist coordinates 
are consistent with zero to within our PN-order accuracy in the q fW metric. This means we can 



employ a similar torus solution method as described in IChakrabartil (|1985l ). A description of 
our modifications to the procedure — including the generalization to our c/>-averaged spacetime — is 
provided in Appendix [A] Note that we now ensure that the equilibriu m solution is found iter atively 
to greater precision, instead of the approximate method described in iDe Villiers et al.1 (120031) which 
has b een used in prior work of the authors in single black hole disk e yolutions (INpble et al.ll2009l . 



2010 ) and by others studying the hydrodynamic circumbinary case (jFarris et al.ll201ll ). We find 



that our procedure produces initial tori that are much closer to equilibrium than the approximate 
scheme. Please see Appendix [A] for more details. 



Previous studies have shown that a gap develops near 2.5a for equal mass binaries (jMacFadyen fc Milosavljevic 



2008: IShi et al.ll201ll ). We aim to study how this gap develops, so we choose to start material out- 



side this radius. We therefore set up a disk with inner edge located at r in = 3ao and pressure 
maximum located at r p = 5ao; from prior experience, r p is approximately the radius at which 
the disk transitions from accreting to decreting since matter must shed its angular momentum to 
fluid elements further out in order to accrete. These outer elements gain angular momentum and 
form a time-averaged decretion flow away from the central potential. We will therefore focus on 
r < r p = 5ao in our analyses. The initial disk extends to r out — 12ao, is isentropic with p/p r = 0.01, 
and is tuned to have an aspect ratio of H/r = 0.1 at r = r p , where H is the density scale height 
defined as the first moment of the rest-mass density with respect to distance from the midplane: 



H 



(p^\e-Tr/2\) 
(p) 



and where the (X) denotes the average over origin-centered spheres: 



(X) 



I 



(18) 



(19) 
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More information about the initial torus and its solution method is given in Appendix[A] In the disk, 
we add random, cell-scale noise to the the internal energy, u, in order to hasten the development 
of turbulence; the random noise is evenly distributed over the range ±5 x 10~ 3 . 

Once the torus is in place on the grid, a surrounding nonmagnetized atmosphere is added 
as our numerical scheme requires us to maintain positive values of p and p. The atmosphere is 
initially static, u l = 0, and in approximate pressure equilibrium: p atm = 1 x 10~ 7 p max (r/M)~ 3 / 2 , 
u atm = 3.3 x l(r 6 u max (r/M)" 5/2 , where 

Pmax and "Umax are — respectively — the initial maxima of 
p and u. We note that when either p or u are found to go below, respectively, p atm or u atm , they 
are set to those atmosphere values without any modification to the magnetic field or fluid velocity; 
this happens very rarely once the disk's turbulence saturates. 

The magnetic field is initialized as a set of dipolar loops that follow density contours in the 
disk's interior. We set the azimuthal component of the vector potential and differentiate it to yield 
B^(t = 0) = in our configuration. The vector potential component is 



= max 



P ~ ^Prnax j ■ 



(20) 



The magnitude of the field, Aq, is set such that the ratio of the disk's total internal energy to its 
total magnetic energy is 100. 



3.3. Grid, Boundary Conditions, and Parameters 

The domain on which the MHD EOM are solved is a uniformly discretized space of spatial 
coordinates {x^} that are isomorphic to spherical coordinates {r, 9,<f>}' 

r(x (1) ) = Me xW , (21) 



0(* (2) ) = \ 



! + (!-£) (2*( 2 ) - l) + U ~ ^) (2* (2) - 1 



(22) 



and cj) = x^ . We set n = 9, £ = 0.87, and 9 C = 0.2. The logarithmic radial coordinates are such 
that the radial cell extents are smaller at smaller radii in order to resolve smaller scale features of the 
accretion flow there. The f-> 8 mapping concentrates more cells near the plane of the disk and 
the binary's orbit, the equator of our coordinate system. Let each grid cell in our numerical domain 
be labeled by three spatial indices that each cover 

[0, ATM _ i], where {N^} are the number of 
cell divisions along each dimension. A cell with indices (i,j,k) is located at (^x f~\ , , where 
Xj = x^ + (j + |) Ax^ n \ The grid we used is completely specified by: x^jp = ln(r m j n /M), 
Ax« = ln(r max /r min )/NW, r min = 15M, r max = 260M, iV« = 300, xf ] = 0, Ax^ = l/N^, 
iV( 2 ) = 160, xf = 0, Ax® = 2vr/iV( 3 ), = 400. 

We chose our resolution and grid extent based upon a number of criteria. First , our 6 and 



resolutions were set in order to adequately resolve the MRI based on guidelines of lHawlev et al 
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(|201 ll ) and ISorathia et al.l (|201ll ) . The radial resolution was chosen to resolve the spiral density 
waves — generated by the binary's time-varying tidal field — by several r adial zones. We find that 
our grid adequately resolves the MRI, as measured by the criteria of lHawlev et al.l (|201ll ) and 
Sorathia et al.l (|201ll ). throughout the domain of interest (i.e. r < 5ao) Vi. Please see Appendix IB1 
for a quantitative description of these resolution criteria and for a demonstration of how well we 
resolve the MRI. 



The radial extent of the grid was inspired by IShi et al.l (|201ll ) and the limits of our near-zone 
PN metric. Since the near-zone P N metric that we use is only valid at distances more than 10M 7 - 



from the black hole with mass Mj (jYunes fc TichsfeOOa ; lYunes et al.ll2006l : IJohnson-McDaniel et al 



20091 ). then — in the equal mass cas e considered here — we have r m ; n > 10f + max(a(*))/2 



5M + ao/2 = 15M. IShi et al.l (120111 ) found that an inner radial boundary located at r m \ n <; 1.1a 
was sufficiently far away from the gap and deep in the potential as to not significantly alter the 
development and evolution of the surface density peak at the edge of the gap. These two constraints 
justify our choice of r m i n = 15M = 0.75ao and suggest that our inner boundary condition may 
begin affecting the gap's evolution when a(t) <^ r m j n /l.l ~ 13. 6M ~ 0.68ao, which occurs after 
approximately t = 51235M in Runln . We set r max = 260M = 13ao to encompass the initial torus. 



All cells are advanced in time with the same time increment (Ax° = At), which itself changes 
in time; Ax° is set to 0.45Ai m i n , where Ai m ; n is the shortest cell crossing time of any MHD wave 
over the entire domain. 



Boundary conditions were imposed through assignment of primitive variables and B l in ghost 



zones. Outflow boundary conditions are imposed at r 



and r 



which amounts to 



extrapolating the primitive variables at th -order into the ghost zones. Additionally, u r is set 
to zero — and u l recalculated — whenever it points into the domain at r = r max . We note that we 
attempted to implement a similar condition on u r at r = r m j n , but foun d it to be unstable during the 
earlies t par t of the simulatio n. Even though it was successfully used in iMacFadyen &: Milosavljevic 
(|2008l ) and IShi et al .1 (|201ll ) , we found that this condition was inconsistent with the tendency of 
negative radial pressure gradients developing ahead of each black hole as it moved around its orbit. 
This pressure gradient moves small amounts of material onto the grid, elevating the density just 
above the floor ahead of the black holes. Even without the special condition on u r at r m j n and 
just using th -order extrapolation of the primitive variables there, we observe insignificant amounts 
(<S 1% of the total) of positive mass flux into the domain there and a nearly flat M(r,t) profile 
over r m ; n < r < 2a with no noticeable artifacts near r m j n (e.g., Figure [7J). 



3.4. Thermodynamics 

Depending on internal properties of the gas (e.g., density, accretion rate), the disk may or 
may not be optically thin, geometrically thin, or have a constant aspect ratio H/r. As these disk 
characteristics are sensitive to the assumed initial conditions and thermodynamics of the system, 
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we therefore must select which kind of disk to model a priori and verify the consistency of our 
assumptions a posteriori. We chose to model a disk with intermediate thickness {H/r = 0.1) in 
order to both address the fact that binary's torque will likely heat the gas efficiently and our 
expectation that the disk will be dense, optically thick and radiating efficiently 

We chose the ideal-gas 'T-law" equation of state to close the MHD EOM: P = (T — 1) pe. 
We set r = 5/3, which reasonably well describes the behavior of a plasma whose specific thermal 
energy is smaller than an electron's rest-mass energy (i.e. it is not relativistically hot). The gas is 
cooled to a target entropy, the initial entropy of the disk, at a rate equal to 



pe {AS | 

Tcool V So 



AS 



Sn 



(23) 



where AS = S — Sq and T coo i = 2ir(r/M) 3 ^ 2 is the cooling time, which we approximate as the 
Newton ian period o f a cir cular equatorial o r bit a t radius r. Our procedure is similar to those 
used by Noble et al. (j2009l ) and Penna et al. ( 2010 ). The term in the parentheses acts as a switch 
ensuring that C c > always, and is zero when the local entropy, S = p/p r , is below the target 
entropy, So = 0.01, which is the constant value used in the initial data's torus. Hence, the cooling 
function should release any heat generated through dissipation since the initial state. We do not 
cool unbound material — i.e. fluid elements that satisfy {ph + 2p m ) ut < —p — since we do not want 
to include cooling that results from application of density or pressure floors. Since C c is the cooling 
rate in the local fluid frame, its implementation in the EOM must be expressed in the coordinate 
frame: 



P u. — £r.Ui 



(24) 



Another advantage of the cooling function is that it provides us with a proxy for bolometric 
emissivity that is consistent with the disk's thermodynamics — unlike a posteriori estimates of syn- 
chrotron and/or bremsstrahlung luminosity that have typica lly been made in numerical relativity 



simulations (e.g.. iBode et al.l (|2010l ); lFarris et all (|2010l . l201ll )). We will use C c to make predictions 



of the total luminosity from circumbinary disks. These predictions are made by integrating C c over 
the domain in the coordinate frame; we expect to verify their accuracy using full GR ray-tracing 
in future work. 



4. Results 

4.1. Approximate Steady State 

At the beginning of both simulations, orbital shear transforms part of the radial component of 
the magnetic field to toroidal, creating a laminar Maxwell stress. Meanwhile, in the same region, 
the magnetorotational instability grows, its amplitude exponentially growing on the local dynamical 
timescale, ~ 500M at the initial inner edge of the disk, r = 60M. The turbulence in the inner disk 



-19- 



reaches nonlinear saturation at t ~ 10000M. Under the combined influence of the initial laminar 
and later turbulent Maxwell stress, matter flows inward (see Figured]). 




100 100 

r [M] r [M] 



Fig. 4. — Color contours of log S(r) as a function of time. The scale is shown in the color bar. The 
black dashed curve shows 2a(t). (Left) Runln . (Right) RunSS . 



4- 1.1. Surface density 

Soon after t ~ 10000M, the inward flow begins to pile up at r ~ 50M, between two and three 
times the binary separation (the dashed line in both panels of Figure [5] marks the location of 2a(t) 
in order to guide the eye). We define the surface density S as 




when we quote it as S(r), that denotes an azimuthal average of equation f)25[) . In later discussion, 
we will sometimes normalize the surface density to So, the maximum surface density in the initial 
condition; in code-units So = 0.0956. In RunSS , S(r ~ 2a) grows steadily for the duration of 
the simulation, but after t ~ 20000M, the logarithmic rate of growth (i.e., dlnS(r)/dt) gradually 
becomes slower and slower. Because a number of azimuthally-averaged properties like S(r) all 
become steadier after t = 40000M, we call the period from then until the end of RunSS the "quasi- 
steady epoch" . For the same reason, we began the binary orbital evolution of Runln at that time. 

Once this quasi-steady state is reached, S(r) rises sharply from the inner boundary at r = 16M 
to r ~ 50M, initially oc r 2,5 , but at late times in RunSS , oc exp(3r/a) (Figure [5]). At first, the 
azimuthally-averaged surface density profile forms a relatively flat plateau at radii greater than 
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Fig. 5. — S(r/a) every 1000M in time from t = 30000M to the end of the simulation. Time 
increases from violet color to red. The dotted curve shows the initial condition. The dashed curve 
shows the average of the colored curves. (Left) Runln , where the time span extends to 53000 M. 
Note that in this simulation a decreases after t = 40000M, so that a fixed value oir/a corresponds 
to a progressively smaller radial coordinate after that time. (Right) RunSS , where the time span 
extends to 76000M. The binary separation is fixed throughout this simulation. 

50M ~ 2.5a, but by t = 30000Af , a distinct local maximum appears at r ~ 50M and persists for 
the remainder of the simulation. This maximum is noticeably asymmetric in the sense that \cK/dr\ 
is always considerably smaller in the disk body (i.e., r > 2.5a) than in the gap region inside r = 2a 
(Figu re El) . This behavior resembl e s clos el y what ha s prey iously been seen in the Newtonian regime 



By construction, the behavior of Runln is identical to that of RunSS up to t = 40000M, when 
the binary inspiral was begun. In fact, at large radius, the behavior of the surface density profile in 
Runln continues to be very similar to that of RunSS even after the binary begins to shrink. Near 
the surface density peak and at smaller radii, however, things change. In Runln , the location 
of the peak moves inward as the binary becomes smaller, and the slope of the disk's inner edge 
becomes noticeably shallower as the inspiral accelerates. Comparing the curve of the dashed line in 
the Runln panel of Figure H] to the curving edge of the colors denoting higher surface density, one 
can see that the location of the disk's inner edge follows the evolution of the binary until shortly 
before the end of the simulation. 

However, speaking in terms of azimuthally-averaged surface density obscures an important 
aspect of circumbinary disks: near and inside their inner edges, their structure is generically far 
from axisymmetric. In Figure El we show S(r, </>) at t = 40000M. As mentioned previously, at radii 
smaller than ~ 2a ~ 40M, there is relatively little matter. The reason this gap forms is that, unlike 
a time-steady, axisymmetric potential, the time-dependent quadrupolar potential of the binary does 
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Fig. 6. — Color contours of surface density in units of So as a function of radius and azimuthal angle 
in RunSS at four different times in two different scales: (Left) Logarithmic color scale emphasizing 
the streams from the disk toward the binary members. (Right) Linear color scale emphasizing 
the growth of asymmetry in the inner disk. In both panels, the times shown are t = 40000A/ 
(upper-left), t = 51963M (upper-right), t = 63926M (lower-left), and t = 75890Af (lower-right). 



not conserve either the energy or the angular momentum of test-particles. Consequently, closed 
orbits do not exist, and torques driven by the binary can rapidly expel s ome matter to the outside, 
while matter on other trajectories can be forced inward (|Shi et al.ll201ll ). As a result, even though 
the rate at which matter enters the gap is comparable to the outer-disk accretion rate, at any given 
time, relatively little matter can be found in the region within ~ 2a, and the matter that is present 
follows trajectories with little resemblance to stationary circular orbits. Instead, a pair of streams 
leave the inner edge of the disk and curve inward toward each member of the binary. Part of their 
flow gains enough angular momentum to return to the disk, but part crosses the inner simulation 
boundary, traveling toward the domain of the binary. 

In add ition, seve r al ten s of binary orbits after matter begins to pile up at r ~ 2.5a, a distinct 
"lump" (as IShi et alJ (|201ll ) called it) forms in the region of the surface density peak. The density 
contrast between this lump and adjacent regions grows steadily in time. Thus, despite the relatively 
slow variation of azimuthally-averaged disk properties during the period we call "quasi-steady" , the 
"lump " continues to e volve secularly. Although RunSS was not continued long enough to see this 
effect, IShi et al.l (|201ll ) found that the eccentricity of the lump's orbit also grows slowly. 
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4-1.2. Accretion rate, internal stresses, and angular momentum budget 




Fig. 7. — Time-averaged accretion rate during four equally spaced segments from t = 30000A/ 
(black) till the end of RunSS . The curves becomes progressively lighter in shade as time advances. 

It is also useful to characterize the global dynamics of circumbinary disks in terms of the radial 
dependence of the net mass flow, i.e., the accretion rate as a function of radius. We show in Figure [7] 
how this quantity slowly evolved during the quasi-steady epoch of RunSS by dividing the time from 
30000M until the end of the simulation at 76000M into four segments and averaging over each one 
separately. The accretion rate is constant as a function of radius only inside the gap region, at most 
times increasing gradually outside r ~ 2a. During the first part of this period, the accretion rate 
rises steadily to radii beyond 5a, but after t ~ 50000M, the accretion rate in the outer disk gradually 
falls. At the end of the simulation, M(r) is actually about a factor of 2 greater at r ~ 3a than 
anywhere else. Averaging over the entire quasi-steady epoch, the rate at which mass passes through 
the inner boundary is a bit les s than half th e accretion rate at r = 5a. Although the first analytic 
theories of circumbinary disks (|Pringldll99ll ) assumed that no accreti on would pass the inner edge 



20081 ) and MHD (|Shi et al. 



2011 



of su ch a disk, N ewtonian simu lations, both purely hydrodynamic (jMacFadven &: Milosavlievic 

, have generally seen leakage fractions of a few tens of percent; 



our fraction is thus only somewhat greater than previously found. 

As mentioned earlier, Maxwell stresses due to correlations induced in MHD turbulence by 
orbital shear dominate angular momentum transport within acc retion disks. Becau se the ratio of 
Maxwell stress to magnetic pressure, 2(BW%)/(B 2 }, is fixed (jHawlev et alJkoilh at ~ 0.3-0.4 



in a point-mass potential (here the notation X^*) denotes the magnitude of the /^-component of 
four- vector X projected into the fluid frame), the stress is linearly proportional to the magnetic 
pressure. A useful measure of the strength of magnetic effects is therefore the plasma (3 = (p) /(B 2 ). 
In most previous accretion disk simulations, this quantity is ~ 100 in the midplane and drops to 
~ 0(1) a few scale- heights out of the plane. We show its dependence on position in the poloidal 
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Fig. 8. — Logio °f the azimuthally-averaged plasma f3 parameter at four times during RunSS 



plane at several times during the quasi-steady epoch of RunSS in Figure to be more precise, 
we show the ratio of the time- and azimuthally-averaged gas pressure to the similarly averaged 
magnetic pressure. As that figure illustrates, the level of magnetization is rather larger than usual 
(i.e., /3 is smaller than usual), but gradually diminishes over time. At t = 30000M, (3 ~ 1 in the 
midplane at r ~ 3-5a and ~ 3 in the region of the surface density peak (r ~ 2-3a); by the end of 
the simulation, it is ^ 10 in the disk body for the whole range 2a < r < 5a and reaches as much as 
~ 30 in the lump. 



Ever since the work of lShakura k, Sunvaevl (|1973l ). it has been popular to measure the vertically- 
integrated, azimuthally-averaged, and time-averaged internal disk stress in units of the similarly 
integrated and averaged pressure. In order to avoid unphysical pressures found in the unbound 
regions, we computed the stresses and pressures only in bound material. Outside r ~ 4a, where 
the disk resembles an ordinary accretion disk, we find that the Maxwell stress alone has magnitude 
~ 0.3-0.5 in these units. This is roughly 3-5 times larger than the stress levels found in general 
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relativistic simulations of MHD flows in the Kerr metric (jKrolik et al.ll2005l ). In the gap region, the 
ratio of Maxwell stress to pressure r ises about a fac tor of 2, while the Reynolds stress in the gap 
rises dramatically (as also found by IShi et al.l (|201ll )). These large Reynolds stresses are entirely 
due to the strong binary torque, which pushes part of the inflowing streams back out to the disk 
with additional angular momentum. 

An overview of angular momentum flow in the system can be gleaned from Figure [9j in which 
we show the radial derivatives of the time-averaged angular momentum fluxes integrated on shells, 
i.e., the time-averaged torques due to the several mechanisms acting. Several important points 
stand out in this figure. The first is that the binary torques are delivered primarily in the gap 
region a <J r <^ 2a. The torque density dT/dr peaks at r ~ 1.45a, and the region surrounding that 
peak dominates the integral over all radii. Moreover, all these torques are positive in net, but they 
are locally negative both at small radii (r <^ a) and at large (r ^ 1.9a). Thus, most of the angular 
momentum the binary gives the disk is delivered in the gap, where the gas density is yery m uch lower 



than in the disk proper. This point has previously been emphasized by lShi et al.l (|201ll ). Second 



that angular momentum is conveyed to the disk proper by fluid flows, i.e., Reynolds stresses. That 
is why the Reynolds stress is large and positive from r ~ 1.8a to r ~ 2.5a. Outside those regions, 
Maxwell stress, which always acts so as to remove angular momentum from the gas and carry it 
outward, dominates the internal stresses. Finally, the net angular momentum change at any given 
radius is generally positive in the inner disk because matter continues to pile up between r ~ 2a 
and r ~ 5a throughout the simulation. 



4-1.3. Disk thickness 

We close this section by commenting on the disk thickness H/r [defined in Equation Q18|) ]. a 
parameter that will play an important role during the period when the binary orbit evolves. Our 
initial data and cooling function were chosen so as to keep H roughly constant over time at a fixed 
ratio to the local radius: H/r ~ 0.1. However, although the gas temperature stayed very close to 
the target entropy at all radii r > 2a, and the ratio H/r did stay nearly independent of radius, 
its value first rose to ~ 0.15 and then fell slightly (to ~ 0.12 by the end of the simulation). The 
departure from the prediction of simple hydrostatic equilibrium was proportional to how much the 
magnetic pressure contributed to support against the vertical component of gravity. 



4.2. Binary Separation Evolution 

At t = 40000M in Runln , we began to evolve the binary orbit, letting it compress as gravita- 
tional radiation removes its orbital energy. The rate of orbital evolution is extremely sensitive to 
separation: a/a oc a -4 when a/r g S> 1. Consequently, even at the relatively small initial separation 
assumed here (oq = 20M), orbital evolution is comparatively slow at first. However, it acceler- 
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Fig. 9. — Radial derivatives of the angular momentum flux due to shell-integrated Maxwell stress 
in the coordinate frame (red), the angular momentum flux due to shell-integrated Reynolds stress 
in the coordinate frame (green), and advected angular momentum (gold). Torque densities per unit 
radius due to the actual binary potential and radiation losses are shown by blue and cyan curves, 
respectively. The net rate of change of angular momentum d r dtJ (solid black). All quantities are 
time-averaged over the quasi-steady epoch in RunSS . 

ates dramatically after t ~ 50000M. By the end of Runln (t = 54000M), a /a is quite rapid, and 
a ~ 8M, small enough to make our PN expansion problematic. 

While the binary orbit changes relatively slowly, the inner edge of the disk moves inward in 
pace with the change in the binary separation, staying close to ~ 2a{t) (as shown in Figures H] 
and [TO]) until t ~ 50000M. However, as the orbital evolution becomes more rapid (after t ~ 
50000M), although the inner edge of the disk continues to move inward in terms of absolute 
distance (Figured]), it begins to recede in terms of r/a(t) (Figure HO]) . At the end of the simulation, 
the disk edge has moved in to ~ 20M, but that is ~ 2.5a. Simultaneous with this evolution, the 
slope of the inner edge also becomes gentler (Figure [5J. In other words, the contrast between the 
surface density in the disk body and in the gap weakens, particularly when considering the outer 
part of the gap. As shown by the RunSS panel in Figure \10\ none of this adjustment (in r/a(t) 
terms) occurs without binary evolution. 

Another view of this process may be seen in Figure [TTJ In that figure, we see the way matter 
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Fig. 10. — Color contours of log £(r/a(i)). The scale is shown in the color bar. (Left) Runln . 
(Right) RunSS . 

accumulates in the inner disk over time, at first during the quasi-steady epoch and later during the 
binary orbital evolution of Runln . The left-hand panel shows what happens when referred to an 
absolute radius scale. When the binary begins to shrink, the quantity of matter found at small radii 
grows abruptly, particularly in the original gap region: the amount of mass inside r = 40M almost 
doubles, and the mass inside r = 20M increases by a factor of 5 during the period of binary orbital 
evolution. The right-hand panel shows the same events from a different point of view. In this 
figure, we see that the mass enclosed within small multiples of a(t) declines rapidly as the binary's 
shrinkage accelerates. For larger multiples (e.g., 3a and 4a), the mass enclosed continues to rise for 
a while after binary orbital evolution, but eventually drops once the compression becomes rapid. 
In particular, the mass within the gap region (i.e., r < 2a(t)) falls by roughly a factor of 40 during 
the period of orbital evolution, although this ratio is in fact a bit ill-defined because 2a(t) is almost 
at the simulation's inner boundary by the end of the simulation. 

The accretion rate behaves differently. It falls (see Figure fl~2|) from ~ 30000M-40000M, even 
before the binary begins to compress. Without binary orbital evolution (RunSS ), it levels out from 
~ 40000M-50000M, before declining more gradually from ~ 50000M until the end of RunSS at 
~ 76000M. In Runln , the onset of binary evolution at t = 40000M leads to a continuing decrease 
in the rate at which mass flows through the inner boundary that levels out only after ~ 50000 M. 
Although the accretion rates in the simulations with and without binary orbital evolution decline at 
different times and at different rates, the final accretion rate in Runln , when the binary separation 
has shrunk to 8M, is only 20-30% less than at the same time in RunSS . 
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Fig. 11. — Mass enclosed within several sample radii as functions of time in Runln . From bottom 
to top, the radii are r/a = 1, 1.5, 2, 3, and 4. (Left) For fixed a = a^. (Right) For time-dependent 
a(t). 

Another consequence of the changing relationship between disk material and the binary is a 
diminution in the integrated torque when the binary compresses (Figure [T3|) . During the initial 
slow stages of energy loss due to gravitational wave emission, the binary continues to exert nearly 
as much torque on the disk as in RunSS , in which the binary orbit does not change at all. However, 
once the orbital shrinkage begins to accelerate, the torque plummets; at the end of Runln , it has 
fallen to ~ 1/5 of the value at that time in RunSS . The greater part of this diminution in torque 
is due to the fact that at this stage in the binary's evolution, its separation diminishes so rapidly 
that the region between a and 2a, where most of the torque is expressed, moves inward faster 
than the matter can follow. There is consequently much less matter on which these torques can be 
exerted. The connection between available matter and torque is shown clearly in the right-hand 
panel of Figure [T3"l in which one can easily see that for nearly the entire inspiral the torque density 
at the location of its maximum (r = 1.45a(i)) is almost exactly proportional to the surface density 
there. However, there is also a smaller part due to an artifact of the simulation. Its inner boundary 
lies at r m ; n = 15M. As soon as a(t) becomes smaller than 15M, part of the region in which the 
binary torque is applied is no longer in the problem volume, so we cannot calculate any torque 
occurring there. As shown in the right-hand panel of Figure \13\ this effect becomes significant at 
t ~ 5.2 x 10 4 M, when a(t) ~ 13M. By the end of Runln , a ~ 8M, so that nearly the entire region 
where the torque is exerted (a ^ r <^ 2a) has left the problem volume. At that point, even if there 
were significant matter there, our calculation can neither say what its mass is nor what torque it 
feels. 
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Fig. 12. — Accretion rate through the inner boundary of the simulation as a function of time. (Left) 
Runln . (Right) RunSS . 



4.3. EM Luminosity: Magnitude, Modulation 



We define the (coordinate frame) cooling rate per unit radius of the disk by 

/— g dO d(f)C c ut. 



dL 

dr 



(26) 



During the approximate stationary state, it is best described in terms of two separate regimes. As 
shown in Figure HU at large radius (r ^ 2a), it is very well described by a power-law, dL/d(r/ao) ~ 
5 x 10 _4 (r/ao)~ 2 Xo a o- At around r ~ 2a, the cooling rate per unit radius reaches a local maximum 
and declines inward. This distinction neatly corresponds to two different mechanisms for generating 
the requisite heat: the dissipation of MHD turbulence associated with mass accretion (at large 
radius) and the dissipation of fluid kinetic energy given to the relatively small amount of gas in 
the gap by the binary torques (at small radius). In fact, this identification is confirmed semi- 
quantitatively. In time-steady accretion, the luminosity per unit radius is (3/2)Mc 4 /[(r/r 9 ) 2 GM] 
at radii where the local orbital angular momentum per unit mass is large compared to the net 
angular momentum flux per unit mass. Our disk is never in inflow equilibrium, and this expression 
is not exact when M is a function of radius. Nonetheless, taking it as an estimator, it predicts 



dL 

dr I ao 



4 x 10" 4 (M/0.01)(r/a )~ 2 S ao- 



(27) 



As Figure [7] shows, the mean accretion rate in code units at r = 2a in RunSS was ~ 0.01, while M 
at larger radii is typically similar or perhaps a factor of two greater. Thus, this prediction of the 
luminosity profile on the basis of the time-averaged accretion rate and expectations derived from 
time-steady accretion onto a solitary mass quite accurately matches the actual luminosity profile 
seen in the simulation. 
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Fig. 13. — (Left) Integrated torque as a function of time in RunSS (black) and Runln (grey). Total 
torque is shown by the solid curves; the dotted curve shows torque in RunSS including only the 
radial range r > r m i n ao/a2(t), where 02 (£) is the orbital separation as a function of time in Runln . 
(Right) Surface density (solid) and torque density at its peak, i.e., at r = 1.45a(t) (dashed) in 
RunSS (black) and Runln (gray). 



Integrated over radius, the total luminosity reaches a peak L ~ 5.5 x 10 3 at t ~ 33000M 
(Figure fT5|) . where L is the integrated luminosity in units of GMT,qc. After reaching this peak, L 
falls slowly, reaching ~ 3 x 1CT 3 at t ~ 76000M in RunSS ; averaged over the entire quasi-steady 
period in this simulation, it is 3.8 x 10 -3 . 

The light output from Runln remains very close to that in RunSS until the binary orbital 
evolution becomes rapid at t ~ 50000M. After that time, it falls more sharply, so that by the time 
at which Runln stops, L ~ 2.7 x 10 ; this is, however, still 2/3 the luminosity in RunSS at the 
same time. As the binary shrinks, the radial distribution of the luminosity changes in parallel, with 
the peak in surface brightness moving inward. We attribute the gradual decline in luminosity to the 
gradual decline in accretion rate. The sharp drop in the final stages of binary or bital shrinkage is 



due to the interaction of a boundary effect with genuine dynamics. As shown bv lShi et al.l (|201ll ). 
gas streams flow inward from the inner edge of a quasi-steady circumbinary disk to radii ~ 1.2a, 
where they can be strongly torqued and some of their material flung back outward toward the disk. 
The outward-moving matter shocks against the disk proper at a radius near that of the surface 
density peak, and the heat dissipated in these shocks contributes significantly to the luminosity. 
When the binary shrinks, this mechanism is weakened for two reasons. The inner boundary of our 
simulation (r = 0.8ao) eventually becomes larger than 1.2a(t); when it does, matter is no longer 
thrown outward by binary torques. At the same time, however, it is possible that the retreat of 
the disk's inner edge when measured in terms of a(t) might also lead to weaker inward streams. 



The fact that the energy deposited by binary torques is ultimately radiated in the disk proper 
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Fig. 14. — Luminosity per unit radius averaged over the quasi-steady epoch in RunSS . The dashed 
line shows a logarithmic slope of -2. 

leads to a method of estimating the relative contributions to the total luminosity coming from 
accretion and binary torques. For that reason, and also because the accretion rate diminishes as the 
region of the surface density peak is approached from larger radius, it is a reasonable approximation 
to suppose that most of the luminosity from the region of the surface density peak inward has its 
source in the binary torques. We can therefore estimate the work done by the torques by bounding 
it between L(r < 2ao) and L(r < 3ao). On this basis, accretion would account for ~ 1/2-3/4 of 
the total (i.e., L ~ 1.8-2.9 x 1(T 3 ) and the binary torque for ~ 1/4-1/2 (L ~ 0.9-2 x 10~ 3 ). 

The rest-mass efficiency of this luminosity is comparable to the rest-mass efficiency due to 
accretion that goes all the way to the black hole. Measured in terms of the time-dependent lumi- 
nosity relative to the time-averaged accretion rate through the inner boundary, the efficiency in 
RunSS falls from a peak ~ 0.06 achieved for 20000M £ t £ 45000M to ~ 0.03 at the end of this 
simulation. There are several reasons that this efficiency is so great even though the potential at 
r = 50M is an order of magnitude shallower than the potential at the innermost stable circular 
orbit (the "ISCO"). One is that the accretion rate in the circumbinary disk is roughly twice the 
accretion rate through the inner boundary, so the local accretion dissipation in the disk is boosted 
by that same factor of two relative to the rate at which mass passes the inner boundary. Another is 
that in a conventional disk around a single black hole the dissipation rate in the region just outside 
the ISCO is depressed relative to larger radii because some of the potential energy released is trans- 
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Fig. 15. — Luminosity as a function of time. (Grey) Runln . (Black) RunSS . We note that the 
vertical axis' range does not include zero in order to accentuate the curves' fluctuations. 

ported outward by the inter-ring stresses. In the Novikov-Thorne model (in which the stresses are 
assumed to vanish at the ISCO), almost 40% of the total luminosity is released outside r = 40M 
when the black hole has no spin. This fraction is smaller when the spin is g reater, and may b e 



further reduced to the degree that the net angular momentum flux is smaller (jKrolik et al.ll2005l ) 



Lastly, of course, additional energy is deposited in the disk by the work done by the binary torques. 
Translating the peak cooling rate into physical units gives 

L disk ~ 2.4 x 10 40 (L/10- 3 )M 6 r erg/s. (28) 

Here To is the Thomson optical depth through a disk of surface density So and L is the luminosity 
in code units, i.e., 3-5 x 1CT 3 . In Eddington units, this becomes L&i^/Le — 1.7 x 10 (L/10 _3 )to. 
Thus, for such a system to be readily observable at cosmological distances, it will be necessary 
both for the disk to be optically thick to Thomson scattering and for the mass of the binary to be 
relatively large. As a gauge of what might reasonably be expected, we note that in a steady-state 
accretion disk around a solitary black hole, the optical depth of the disk at r/r g = 20 would be 
~ 2 x 10 3 (a/0.1)~ 1 (ry/m), where r\ is the usual rest-mass efficiency and rh is the accretion rate in 
Eddington units. With this disk surface density, the luminosity would approach that of a typical 
AGN when Mq is at least ~ 1. 

If this light were radiated thermally, the corresponding effective temperature would be 

T cS ~ 4 x 10 4 (L/10- 3 ) 1/4 M 6 " 1/4 r 1/4 K, (29) 



where we have assumed that the radiating area is 2vr(2a) 2 . Thus, it would emerge primarily in the 
ultraviolet for fiducial values of black hole mass and optical depth. 
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Fig. 16. — Fourier power spectrum of the luminosity radiated during the quasi-steady epoch of 
RunSS . The vertical lines represent: the orbital frequency at the surface density maximum (dashes) 
and the peak in the spectrum (dots). 

The luminosity (assumed to be optically thin) exhibits a noticeable modulation as a function 
of time, with peak-to-trough contrast of ~ 5%. Its Fourier power spectrum shows a strong, sharp 
peak at a frequency 1.47f2bin (see Figure [To]) and a weaker peak at 0.26Slbin- The latter is the orbital 
frequency at the radius of the surface density maximum, ~ 2.4a; because the lump is located at 
this radius, we call this frequency Oiump- The former we identify with the rate at which the lump 
approaches the orbital phase of a member of the binary, 2(Qbin — ^lump) = 1.46f2bm- When the 
lump draws near one of the black holes, a new stream forms, falls inward, and is split into two 
pieces, one of which gains angular momentum, sweeps back out to the disk, and ultimately shocks 
against the disk gas. It is this process, whose frequency is 1.46f2bm, that modulates the light curve. 
If the binary mass ratio were far from unity, we expect that the modulation frequency would fall 

to — f^bin ^lump- 
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5. Discussion 
5.1. Comparison to Newtonian MHD 

In many respects, the behavior we fo und in this pos t-Newtonian regime resembles what was 



previously found in the Newtonian limit (|Shi et al.l l201ll ). There is very good agreement in the 



shapes of their azimuthally-averaged surface density profiles, with any contrasts attributable to 
their somewhat different initial conditions. In both cases, during the quasi-steady epoch the surface 
density at the disk's inner edge rises oc exp(3r/a), reaches a maximum at r ~ 2.5a, and then declines 
to larger radii. 

At early times in both, there is a pair of streams leading from the disk edge to the inner 
boundary, which they typically reach at an orbital phase slightly ahead of the nearest member 
of the binary. Both also develop a strong m = 1 asymmetry (a "lump") in the surface density 
at r ~ 2.5a at late times. This asymmetry ultimately causes, in both the Newtonian and post- 
Newtonian simulations, a single stream in the gap to become dominant. Almost the only contrast 
in this regard is that the orbit of the lump developed a growing eccentricity in the Newtonian case, 
but not in RunSS . 

The level of magnetization is likewise qualitatively similar: the mean plasma (3 in the Newto- 
nian case fell from ~ 1 at ~ 6a to ~ 0.3 at r ~ 2a, while the value (averaged over the quasi-steady 
epoch in RunSS ) in our simulations was ~ 1.5 at r = 6a, grew to ~ 2.5 at the surface density peak, 
and then decreased inward. The magnetic stress-to-pressure ratio a in the disk body follows the 
same pattern of close resemblance. It was ~ 0.3 in the disk body in the Newtonian case, and ~ 0.2 
in RunSS . In the gap, the similarity was more qualitative than quantitative: in both cases, it rose 
steeply into the gap, but reached only ~ 0.7 at r ~ a in the PN simulation, whereas it climbed to 
~ 10 in the Newtonian one. 



Most strikingly , the l uminosity estimated by IShi et al.l (|201ll ) scales extremely well to the 



PN case. IShi et al.1 (120111 ) could not directly compute the luminosity because they assumed an 
isothermal equation of state. However, they argued that the work done by the binary torques would 
be delivered to the disk and ultimately dissipated there into heat. Rewriting in our units their value 
for the rate at which the torques did work on the gas gives a luminosity of 0.018G , MS p c(a/r 9 ) -1 / 2 , 
where S p is the surface density at the maximum; for our separation (a = 20M) and our surface 
density at the maximum (~ 0.55 averaged over the quasi-steady epoch in RunSS ), that becomes 
2.2 x 10 _3 GMSoc. This prediction agrees well with the upper end of our estimated range for the 
binary torque share of the luminosity. 



5.2. Comparison to Analytic Estimates of Binary Runaway 



Milosavljevic Phinneyl (120051 ) predicted that at some point well before the merger, the 
BBH should begin compressing so fast by gravitational radiation that internal stresses within the 
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disk would not allow it to move inward rapidly enough to stay near the binary. At the order of 
magnitude level, this breakaway point would be expected to come when the gravitational radiation 
time 

_5_ / a Y (1 + g) 2 



t = ^ U-M (30) 

g 64 \MJ q v 1 

becomes shorter than the characteristic disk inflow time 

t in = a- l {H/r)- 2 {d\n^/d\nr)- x 9,- 1 = a- 1 (F/r)- 2 (dlnS/(ilnr)- 1 (r/r 9 ) 3 / 2 M. (31) 

In these equations fi is the local disk orbital frequency. The logarithmic derivative of the surface 
density enters because spreading of the inner edge is more rapid when it is especially sharp. 

With a typical estimate of the stress level, a ~ 0.01, the binary separation at which and 
tin would match, and the disk and binary might decouple is 

(tt I \ —4/5 
oil) M ' (32) 

where the fiducial radius at which the inflow time was computed is r* = 2a, and we set q = 1. 
Indeed, it was this sort of estimate that led us to choose the initial conditions for our simulation. 

However, scaling to the actual parameters of our simulation leads to a considerably smaller 
predicted value, 

4/5 



10[(dlnS/dlnr)/6]- 2 / 5 (— ) M, (33) 



which is much closer to what is found in Runln . Thus, the substantially stronger magnetic stresses 
than predicted by usual a-based estimates lead to decoupling at a much smaller binary separation. 
Nonetheless, in the end the inward motion of the disk is limited by angular momentum transport, 
so once the magnitude of those stresses are known, adec can be estimated quite well by this means. 

On the other hand, the meaning of the term "binary runaway" should also be made more 
nuanced. As we have seen, the accretion rate through the simulation inner boundary decreases as 
the binary shrinks, but almost the same decrease in accretion rate occurs when the binary does not 
shrink. Moreover, the continuing advance of the disk during the period of orbital evolution brought 
matter rapidly inward. The amount of matter within 30M rose by about a factor of 4 while the 
binary shrank, so that almost as much matter could be found within that radius as had been within 
40M at the beginning of the binary orbital evolution. Thus, acceleration of binary orbital evolution 
does lead to a state in which the surface density at r < 3a is smaller than would be expected if 
the orbital evolution were slower, and this diminution in the mass close to the binary does lead to 
consequences such as sharply diminished torque (as discussed in Section 14. 2\\ and luminosity (see 
Section T4.3p . On the other hand, neither the torque nor the luminosity falls by as much as an order 
of magnitude because the decoupling of binary and disk matter is not complete. Most notably, 
accretion continues at a rate only tens of percent lower than in the absence of inspiral. 
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Continuing accretion into the binary orbital region has a particularly interesting consequence. 
The material in those streams should, just as happens when the binary orbital evolution is slower, 
be captured into orbit around one or the other of the members of the binary. It will then settle 
into two smaller disks, one around each black hole. In the conditions of our simulation, the inflow 
time in the individual black hole disks might be only slightly shorter than the merger time because 
decoupling occurs when the binary separation is not a great deal larger than the ISCO, even if 
both black holes spin rapidly. If the circumbinary disk were cooler than in our simulation, so that 
a dec 3> M, the inflow time in the smaller disks would be shorter than that in the circumbinary 
disk by a sizable ratio: ~ 3 x 10~ 4 at decoupling if the scale heights in the inner disks and the 
outer disk are the same, and even smaller at later times. In either case, there could be interesting 
hydrodynamic interaction between the two smaller disks as the binary compresses. 



5.3. Different Disk Thermal States 

In terms of an ultimate comparison to observations, a larger question is posed by the disk's 
thermal state. As just shown, the binary separation at decoupling scales as (if/V) -4 / 5 , so the disk's 
internal pressure (H oc c s oc (p/p) 1 ^ 2 ) can influence adec- For reasons of numerical convenience, we 
chose parameters yielding a relatively thick disk. Although the factors controlling the saturation of 
magneto-rotational turbulence are still not well understood, a scaling with disk pressure remains 
plausible. If the effective sound speed of the gas were lower, decoupling might occur at rather larger 
binary separation, well outside the domain of relativistic orbits. 

Several factors can influence the actual equation of state of the disk. In ordinary AGN, local 
heating due to accretion can make the disk radiation-dominated inside r ~ 100M when m ^ 0.3 and 
the central mass M > 1O 6 M jKroWJlQfld ). Larger central masses lead to radiation dominance even 



when m is smaller. When the disk surrounds a binary, the local heating should be similar at radii 
r ^ 2a, as shown by Figure HU the diminished accretion inside ~ 4a is compensated by dissipation 
of the work done by the binary torques and delivered to the disk. In those circumstances, the disk 
scale height for r ^ 2a is independent of radius, giving an aspect ratio H/r ~ (3/2)(m/r/)(r/r 9 ) _1 , 
where m is now the local accretion rate in Eddington units. Thus, our aspect ratio of ~ 0.15 would 
correspond to a nominal rh ~ 0.4(r/40r 9 ); that is, this m(r) is the mass accretion rate at r that 
would, if it reached the black hole in a flow with r\ = 0.1, produce that fraction of an Eddington 
luminosity. 

Given the relatively large leakage fraction through the inner edge of the circumbinary disk, the 
accretion rate onto the two black holes will in general be smaller than the accretion rate in the disk 
by only a factor of a few. They might therefore generate a sizable luminosity with a spectrum not 
too different from that of a generic AGN. Because the density in the gap is considerably smaller 
than in the disk, this luminosity may irradiate the disk, particularly if it is relatively thick. The 
inner edge of the disk could be heated by this means, as well. 
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Thus, there remains considerable uncertainty in the thickness profile of a circumbinary disk 
surrounding a relatively compact binary black hole. The particular thickness we have simulated 
lies within that range of uncertainty, but may not be generic. 



5.4. Distinctive EM Signals 

Given a sufficient external mass supply rate, the luminosity from the circumbinary disk alone 
could be great enough to be detected, even from a cosmological distance (Section I4.3|) . However, 
until the binary separation becomes as small as that considered here, the luminosity from matter 
accreting onto the two individual black holes would dominate the circumbinary luminosity by a 
large ratio. In many respects, a binary black hole with a 3> 10M should strongly resemble a 
conventional AGN. However, when the separation is as small as the ~ 20M scale studied in our 
simulations, contrasts with ordinary AGN continua might occur due both to the additional heating 
near the inner disk and the gap from ~ 2.5a down to ~ a/3 in the range of radii in which a thermal 
disk exists. The supplementary heating due to the binary torques will increase the luminosity of 
the portion of the disk near r ~ 2a, but the temperature in this region is smaller than the hottest 
part of the accretion flow by a ratio ~ (1.5nsco/2a) 3 / 4 , where the factor 1.5 multiplying the ISCO 
radius is meant to account approximately for the displacement of the temperature maximum from 
the ISCO. Consequently, the additional luminosity will appear at rather longer wavelengths than 
the peak of the thermal continuum. On the other hand, radiation from the gap region will be very 
different from what might be expected from a conventional disk in those radii. The dissipation 
rate is much smaller because the motions are laminar, not turbulent; moreover, whatever light 
does issue from that region is unlikely to be effectively thermalized, and would therefore emerge 
at considerably shorter wavelengths. Thus, the luminosity at wavelengths intermediate between 
those characteristic of the innermost part of the disk and those characteristic of r ~ 2a would be 
significantly suppressed. 

The periodic modulation in the heating rate of the circumbinary disk that we have found 
might make its emission easier to isolate. The key question governing that "might" is how effec- 
tively optical depth in the disk blurs the modulation. Our cooling function, which operates at a 
characteristic rate ~ il(r) filters out variations on timescales <C Q , but optical depth would im- 
p ose a rather more severe upper bound on the maximum effective frequency of variation. According 
to lKylafis fc Klimisl 1)19871 ) . a constant-density sphere of radius R and optical depth r suppresses the 
amplitude of a periodic signal of angular frequency oj injected at its center by a factor ~ 3c/(Rtuj) 
when lvRt/c ^> 1. To apply this estimator, we suppose that a stratified disk segment with scale 
height H can be approximated by a homogeneous sphere of radius R = H. For binary separation 
a, the relevant orbital radius is ~ 2.5a and the signal frequency u = 1.46f&bin ~ 1.5(GM/a 3 ) 1 / 2 for 
total binary mass M. As shown in Figure [6l the surface density in the lump grows to be ~ Eo, 
so we take t = tq. The suppression factor is then ~ 0.024(r /1000)- 1 (a/20r 9 ) 1 / 2 (i?/0.15r)- 1 . In 
other words, when the relevant region of the disk is optically thick, the luminosity in the modu- 
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lated component is independent of surface density, so that its fractional modulation decreases as 
the luminosity increases. 

As we have discussed in section [5TB"l there is also considerable uncertainty in H/r, even given a 
disk surface density. If the disk thickness and optical depth were determined by the considerations 
of conventional time-steady accretion flows around single black holes, the characteristic cooling time 
tH/c can be identified with (afi)" 1 . The fluctuation suppression factor could then be estimated by 
~ 2>a£l/u. Here the relevant orbital frequency is 0(2. 4a) = 0.26r2bin> while 0.74f2bin < w < 1.46f2bi n 
(the upper limit applies in the equal-mass case, the lower limit when the masses are very unequal). 
Thus, the suppression factor estimated in this way is ~ 0.1(a/0.2) for equal black hole masses, 
rising to double that in the limit of very unequal masses. However, this estimate is made uncertain 
by the fact that the disk in the vicinity of the surface density peak is certainly not in a state 
of inflow equilibrium. In addition, just as for the other estimate, the association of the periodic 
modulation with the lump means that any estimate based on assumptions of axisymmetry likely 
underestimates the local optical depth. 

Both estimates suggest that the modulation will be suppressed by at least a factor of several, 
but both are also subject to considerable uncertainty, making the actual outcome unclear. It is 
worth pointing out that in the event the modulation is detectable, the period of the modulation 
would allow an estimate of the binary orbital period. When the binary mass ratio is unity, the 
binary orbital frequency is 0.68 times the frequency of the modulation; as the mass ratio departs 
from unity, the binary orbital frequency should rise toward ~ 1.36 times the modulation frequency. 

Finally, we remark that our predictions of EM signals f rom circumbinary disks around merging 



black holes are complementary to those previously made (jBode et al.ll2012l ) in two ways. In the 
previous work, the period of EM emission began when the binary separation shrank to 8M; that is 
when our calculation ends. In addition, that effort expressly excluded the disk proper, which they 
defined as r > 16M; in our work, that is the location of the overwhelming majority of the emission. 
Our effort also differs from previous work in this area in that we explicitly include radiation losses 
in the gas's energy equation and also tie the rate of radiation directly to the instantaneous local 
thermodynamic state of the gas (albeit in only a formal way). In addition, our discussion of the 
observability of periodic modulation in the lightcurve takes into account possible suppression of the 
variation due to optical depth in the source. 



6. Summary 

By describing the binary black hole spacetime at separations of tens of gravitational radii 
through the PN approximation, we have been able to simulate many orbits of fluid motion around 
such a system in fully relativistic MHD. In so doing, we have demonstrated that the qualitative 
properties of circumbinary disks in such a regime are well described by an extrapolation from their 
properties in the Newtonian limit: Matter piles up at ~ 2.5a, while smaller radii are largely cleared 
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of mass; nonetheless, accretion continues through the inner gap, albeit reduced by a factor of a few 
from the rate at which it is supplied at larger radii. 

At the same time, however, we have also investigated the initial stages of strongly relativis- 
tic behavior in the form of the disk's response to binary orbital evolution by gravitational wave 
emission. By carrying the disk's evolution through the transition from the epoch in which its char- 
acteristic inflow time is short compared to the binary evolution timescale all the way to the epoch 
in which the binary evolves much faster than the disk, we have established the time at which the 
binary "runs away" from the disk, and more importantly, the degree to which it does so. This de- 
coupling causes a drop in both the torque the binary exerts on the disk and in the disk luminosity. 
However, a sizable fraction of the accretion rate at large radius continues to makes its way to the 
binary throughout this period. 

The binary separation at which this decoupling occurs is rather smaller than commonly esti- 
mated, largely because the internal stresses produced by MHD turbulence in the disk are consid- 
erably greater than typical applications of the a-model had guessed. The actual value of ad cc is 
sensitive to the disk's thermodynamics to the degree that the absolute level of the internal stresses 
are proportional to the disk's internal pressure. Because accretion continues, luminosity released 
when the accreting gas reaches the black holes may illuminate the disk and heat it. This sort 
of feedback has the potential to keep the disk's inflow rate high, self-consistently sustaining the 
accretion rate. 

Given the sort of accretion rates associated with AGN, the inner regions of circumbinary disks 
around binary black holes with separations of tens of gravitational radii can be almost as bright as 
AGN, although there may be identifying alterations in the shapes of their optical/UV continua. 

We have also shown that the work done on streams passing from the inner edge of the circumbi- 
nary disk through the evacuated gap around the binary is carried back to the disk and dissipated 
there. Because the disk generically develops a non-axisymmetric density distribution at ~ 2.5 bi- 
nary separations, the dissipation rate is modulated periodically with ~ 5% fractional amplitude. 
In the right circumstances, this modulation might be detectable, although optical depth in the disk 
is likely to diminish its fractional amplitude, particularly when the accretion rate is high enough 
to make the system luminous. If this modulation can be detected, its period would provide an 
estimator of the binary's orbital frequency with factor of 2 accuracy. 

A number of our results may be sensitive to the particular parameters chosen, most importantly 
equal masses in the binary, spinless black holes, and perfect alignment between the orientation of 
the binary's orbital angular momentum and the disk's angular momentum. Moreover, these choices 
can interact: for example, black hole spins oblique to the gas orbital plane can induce changes in 
that plane. Future work exploring a variety of choices for these parameters may reveal additional 
effects. 
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A. Hydrostationary Torus Solutions in General Axisymmetric Spacetimes 



Here, we describe a method for calculating axisymmetric non-magnetized gas distributions 
supported by pressure gradients and rotation within an axisymmetric spacetime. We will assume 
a general spacetime with Killing vectors (d/d(j)) a and (d/dt) a that can be expressed in the simple 
form (in coordinates similar to spherical Boyer-Lindquist coordinates): 



which means that the inverse metric is 
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where A = gf^ — gttg<j>4>- We have verified that the ^-average of our PN spacetime, g^, has this 
form to within the accuracy of our PN procedure. 

The initial state of the simulation consists of matter in axisymmetric hydros tatic equilibrium 
with a specific angular m omentum pr o file, £. We start from the discussion of iDe Villiers et al 



(|2003l ). which is based on IChakrabartil (|1985l ) and other citations mentioned therein. The disk is 
centered about the equator of the black hole's spin; we will eventually assume that it is initially 
isentropic. The time-independent and axisymmetric Euler-Lagrange equations reduce, essentially, 
to 

nh i o 
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where the angular frequency — Q = u^/u t — is not a simple function of the specific angular momentum- 
£ = —Ufk/ut. The 4-velocity, u^, in our symmetry has zero components: u r = vP = u r = u$ = 0. 

-1, that 



One can show, from the normalization condition u^u^ 1 
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where r] and q are yet to be determined parameters, and A is defined by 
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We can eliminate Q from this system by combining equations (|A6j) and (1A7I) to yield a non-linear 
algebraic equation for I = £{r, 9) in terms of the metric: 



where 



or 



R(£) = tf* [£ 2 + \ 2 {£)] - g a £ - g^£\ 2 (£) = , 
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= r]\ 2 ~ q , 
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Also, we can show that 
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= v -2/(q-2) e q/(q-2) 
= kf, 
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where k = r] 2 ^ q 2 ^ and ( = q/ (q — 2). 



A 



Typically, one "solves" equation (IA8I) by approximating A with its Schwarzschild value: 

i. For our ^-averaged 
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(De Villiers et al. 2003; Noble et al 



2009; 



2011 



.Farris et al. 

spacetimes, the Schwarzschild approximation is not so gooc|l|. Therefore, we need a better solution. 



7 When using the approximate method, we found the disk to undergo a low frequency breathing mode that 
dominated the early evolution of the disk. 
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We can solve this equation to roundoff precision by using a Newton-Raphson scheme. To do so, we 
will need to know dR(£)/d£: 
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where dX 2 /d£ = 2A 2 /[(2 — q)£]. Let us come back to equation ()A3p . A solution to this equation 
yields our disk solution. The solution process involves integrating it from the inner disk's edge - 
located at r- m — to a point within the disk: 
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dh 

With the boundary condition, h- u 
h = h(r,0): 



d(u t ) 



(«t)" a 

h(r in ,9 = tt/2) 



+ 



l-k£t +1 



d£. 



(A13) 



h 



U t \nf{£h 



1, one can solve this integral equation for 
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where Ut is given by equation (|A4p 
^( r in>7r/2) is a boundary value, and f(£) = 1 — k£^ +1 



ut(r,6)f(£(r,e)) ' 
is found using Newton-Raphson on equation (|A8p . 
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We want a distribution that resembles a torus, which has a pressure maximum at some radius 
and has finite extent. The parameters {£; n , q, rf\ determine whether we get such a solution. We 
would like to replace one of the degrees of freedom with r„, however, there is not a closed-form 
solution for r p in terms of any of the original parameters. We know that the fluid attains the 
Keplerian angular momentum at the pressure maximum (r p ) as the pressure gradient must be zero 
there. It means that £\ n should be super-Keplerian at the inner edge (n n ), and £ out should be 
sub-Keplerian at the outer edge (r out ). We therefore know that £- m > ^(Hn), 4>ut < ^x(^out), and 
£p = £K( r p)i where £k is the Keplerian specific angular momentum (see the next section below). 
However, we only have two free parameters, and now have three constraints (if we specify all -^in, ^out 
and £ p ). It may be possible to change equation (|A9P to look like: 

,2-q 
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and then find Ao with this third constraint. Using these three constraints, however, does not yield 
a closed- form solution for q, rj, and Ao- Therefore, another Newton-Raphson procedure would be 
required. Hence, we relax the constraint on 
Using Aq = 0, we find that 



out < ^ftr(^out)) and let r out be a result of our procedure. 
n log (£ iD /£ p ) 



log(A in /A p ) ' 



Ap 



(A16) 
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where X p = X(£ p ,r p ) and A; n = A(^;, 



given by equation (|A7[) . 



We follow the solution process based on one described in IChakrabartil (|1985l ). and is the 
following: 
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1. Chose values of (ri n ,r p ,£i n ), and derive q and r) using equations (|A16P and (|A17P : 

2. Calculate A (and then £) at a new (r, 8) via Newton-Raphson (using equations (1A8|) and 
(lATOl) ) close to the previous location, so that the old location's value can be used as a seed 
to the Newton iteration to successfully find a solution; 

3. Calculate u t via equation t|A4j) then h at (r, 6) via equation (|A14p : 

The solution process progresses throughout (r, 9) space until the boundary of the disk is found, 
where h(r,9) = 1. The path we take starts at (r,6) = (rju, 7r/2), moves in increasing r along the 
9 = tt/2 line, and we test to see if /i is increasing. If h is not increasing at first, then we stop and 
try a different set of parameters. If h is increasing at first, we then proceed until the outer edge of 
the disk is reached; this radius is denoted as r out . Then, Vr € [ri n ,r out ], we start from the 9 = tt/2 
solution and proceed backward and forward in 9 along constant r to find h(r,0). 

The initial torus solution used herein is parameterized by l- m = 8.743, rj = 1.961, q = 1.642, 
r in = 3a = 60M, r p = 5a = 100M, r out = 11.75a = 235M. 



A.l. General Keplerian Velocity 



The stationary torus solution described in Section [X] requires the Keplerian, or circular equa- 
torial geodesic orbits, of the spacetime. Since we do not calculate g^ u in closed form, we require 
equations for these orbits based on a generalized metric of the form Equation (|A1|) . Here, we state 
the equations governing the Keplerian orbits. 

Keplerian orbits in our spacetime have 4-velocity, u^ 1 = [u*,0, 0, u^} = it*[l, 0, 0, Qjc]- is 



found from the r-component of the geodesic equation, du r jdr 
yields 
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$1jc _ and £Ik + are the prograde angular velocity and retrograde angular velocity, respectively. The 
4-velocity components are found by the normalization condition: 
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The Keplerian specific angular momentum, £jc, is found by the relation between £ and 0: 
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Fig. 17. — Change in total energy relative from its initial value versus time of non-magnetized tori 
from different runs. The first run used the so-called "approximate method" and was evolved in a 
PN BBH spacetime (solid); the second run used our new and more accurate procedure, but used 
the same spacetime as the first (dashes); the last used the same disk from the second run, but was 
evolved in the (/>-average of the other runs' spacetime (dots). The curve for the last run (dots) 
oscillates with amplitude ~ 10 -8 , which is why it appears consistent with zero in this figure. All 
the disks share the same parameters: r ia = 2ao, r p = 4ao, ao = 30M. 



B. Resolution Requirements 

Our grid resolution was chosen to adequately resolve the MRI, to resolve the spiral density 
waves generated by the binary's potential, and to involve cells that are nearly cubical. We discuss 
each choice in turn. 



Many r ecent studies have explored the resolution dependence of glob al MHD accretion disk 



simulations ( Hawlev et al. 



2011 



Sorathia et al 



2011 



Shiokawa et al.1 \20m . lHawlev et al.1 (|201lh 



found that many global properties of the disk nearly asymptote with increasing resolution once the 
following criteria are satisfied: 
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iV (3) ;> 790 



0.1 
H~fr 
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(B2) 



where is the number of cells per scale height, H, Q( z ' > 10 is the recommended quality factor 



of the simulation, f} z = (p) / (\y/g zz B z \ 2 ) . We use spherical coordinates, so is needed instead: 

H H/r 
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(B3) 



where Njj/ r is th e number of cells in the poloidal direction per scale height. We see from prior 



simulations (e.g., iNoble et al.l (j2010j)) that (3 ~ 10 and /3 z //3 ~ 50 are reasonable for a disk in its 
asymptotic steady state, suggesting that Njj/r > 36. The initial condition values of /3 ~ 100 and 
/3 2 //3 ~ 1, however, yield a weaker constraint (N H i r > 16) on the resolution. Thus, we setup a grid 
such that Njj/ r — 36 with H/r = 0.1, our simulation's scale height. This is satisfied by the 
discretization described in Section [3731 We also note that our condition satisfies the recommendation 



of N H/r > 32 by ISorathia et al.l (|201ll ) 



The more severe constraint is on the azimuthal symmetry. Both lHawlev et al.l (|201ll ) and 



Sorathia et al.l (|201ll ) suggest that past simulations under-resolved the azimuthal direction and 
that one should cover the full azimuthal range <f> E [0, 2ir] instead of assuming quarter- or half-circle 
symmetry. Since A(f> limits the time step size, we were only able to afford N^ 3 ' = 400 as anything 
larger was impractical given our computational resources at the ti me. We were optimistic with this 
resolution, however, since the thinnest run of INoble et al.l (|2010i ) failed to satisfy Equation (jB2]) 



yet still resolved the MRI with > 25 throughout most of the disk's body. 

We demonstrate how well Runln and RunSS resolve the MRI in Figures [18] - \19\ where we 
show mass-weighted averages of the Q (2) and Q (3) MRI quality factors: 



Q 



The averages were made over in the following way: 
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A mass-weighting is used to calculate {Q l ) P in order to bias the integral over the turbulent por- 
tion of the disk (the disk's bulk) rather than the laminar regions (e.g., corona, funnel). We find 
that the constraint, i.e. (Q^) p > 10, is satisfied for all times and regions in either Runln or 
RunSS except for the densest parts of the lump at late times in RunSS . Similarly, the con- 
straint, i.e. (Q^)n > 25, is satisfied for all times and regions in either Runln or RunSS except for 
in the lump at late times in RunSS . We further note that (B r B<f ) )/(p m ) ~ 0.3 — 0.35 when averaged 



found in resolution studies about point masses (Blackman et al. 


2008; 


Guan et al. 


2009 


Davis et al. 


2010; 


Simon et al. 


2011; 


Hawlev et al. 


2011; 


Sorathia et al. 


2011: 


Shiokawa et al. 


2012). 
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We also aim to resolve the spiral density waves generated by the binary's time-dependent 
quadrupolar potential. This means that we need about ~ 10 cells per wavelength of the sound wave 
generated by the binary, = 2 ir c s /Qun, where c s = (H/r) r VLk is the speed of sound. We use the 
Newtonian approximates fibin — a~ 3 / 2 and — r~ 3 / 2 to simplify A^: A^ ~ 2ir (H/r) r (a/r) 3 / 2 . 
We want to resolve A^ in <j) and r out to r p , which means that we want to satisfy another quality 
condition: A^/Ar, X^/rAcf) ~ Qd, where Qd will be the target number of cells per spiral density 
wavelength. Using the grid specifications described in Section 13.31 it is easy to show that the 
resolution constraints become: 

We find that the spiral density wave criterion is stricter than the MRI criterion when 

Q,Q 3/2 >/3 1/2 Q (3) . (B8) 
Again, we did not satisfy the constraint on azim uthal reso l ution with our choice of N (3) = 400. In 



this case, we were reassured by evidence found bv lShi et al.l (|201ll ) that found that the spiral density 



waves were extended over a large azimuthal extent and required far fewer cells than expected in 
that direction to resolve. In practice, we find that spiral density waves are short-lived as they 
propagate through the turbulent, shear flow of the disk; they are hardly ever seen in snapshots of 
intrinsic quantities past r ~ 3ao- 

Since gridscale dissipation scales with the ratio of the cell extent to characteristic length scales 
of the physical quantities, cells that are oblate may effectively lead to anisotropic dissipation. 
Because physical dissipation mechanisms are isotropic, this effect could lead to unphysical artifacts. 
For this reason, we attempt to make the cells within the bulk of the disk — where most of the 
dissipation occurs — as isotropic as possible. Our runs use grids with Ar : rA9 : rAcf) :~ 3 : 1 : 3 as 



measured in the 9 = tt/2 plane. ISorathia et al.l (|201ll ) suggested that Ar = rAcft and A<p/ AO < 2; 



we satisfy the former, and violate the latter by a slim margin: A<p should be just 2/3 times the size 
we use. We note that the poloidal extent increases off the equator, so that the cells become more 
cubical at larger \6 — n/2\. 



C. Mass, Energy, and Angular Momentum Budgets 

Space and time gradients of the accretion flow's extensive quantities, mass (M), energy (E) 
and angular momentum (J) are fundamental for understanding the disk's evolution and structure. 
In stationary spacetimes, M, E, and J are all conserved. In our (t, </>)-dependent spacetime, E and 
J are no longer strictly conserved. We describe here how several functions used throughout the 
paper are derived from the evolution equations of mass, energy and angular momentum. 



-46- 



C.l. Angular Momentum 



We begin this section with the angular momentum equation beca use of its import to accretion 
physics. We follow the notation and derivation procedure outlined in lFarris et al.l (|201ll ). 



An extensive quantity, J, is the integral over the spatial volume of the time component of 
the its associated current, j^: J = J j* \J—g dV, where dV is the spatial volume component in 
the spacelike hypersurface (e.g., drdOdfi). We are interested in the azimuthal component of the 
momentum as that is the dominant component of the gas' and the binary's momenta. We therefore 
recognize that j M = T tM u (j) u , and <p u = {d^Y = dx v /dcj) = [0,0,0, 1] in spherical coordinates, which 
is what we use. We wish to calculate d 2 J/dtdr . If J is locally conserved perfectly, V M j M = 0. In 
our case it will not be conserved exactly, and exploring the radial gradient of its volume integral 
will help us understand how MHD stresses and the binary's gravitational torque compete over the 
run of the flow. This quantity is: 



^ f (^.f ) ^drdOd^ = d r f (d„ 
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where the last equality results from the fact that j is zero on the axis, and j^(<p = 0) = 
On the other hand, we know from the stress-energy EOM — V ^T^j, = — T v — that 
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where the torque density, dT/dr, can be expressed as 

We remind the reader that J- v is the radiative cooling flux (see Section [3] for details). 
Therefore, equating the two equations ()Cip and (|C2[) . we have 
drdtJ = § - - d r {T r ^} 

where we have used here the shorthand 



{X} = J ^g-XdOdcj) =(X)JV=9 



de, 



(Cl) 
= 2vr). 



(C2) 



(C3) 



(C4) 



(C5) 



Also, M r ^, R r 'tj,, and A T $ are — respectively — the Maxwell (MHD) stress, Reynolds stress, and 
advected flux of angular momentum. We note that M^ v = 2p m u^u u + p m S^ u — b^b u is the EM 
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part of T^v, while (R^ u + A^ u ) = Th^ v = phu^u u + p5^ v is the hydro dynamic part. The Reynolds 
stress alone is more complicated to calculate as we have to find the perturbation from the mean 
flow: 

R r fj, = ph 5u r Su^ , (C6) 

where 

5u^ =u^ - {put 1 } I {p} . (C7) 

We note that we include the enthalpy as it technically contributes to the stress; its contribution 
is insignificant, however, for our relatively cool flow. The quantities {R^ u } and {A^ u } are not 
calculated during the simulation, but found approximately from other shell-integrated quantities 
we do calculate; {M^ l u }, {J 7 ^}, and dT/dr are calculated as stated above during the run. One can 
easily show from equations (|C6|) and (|C7[) that 



{R r ,} = {ph5u r 5 U( p} ~ {T H r 4-{^ r 4 (C8) 
where Tjj r \ is the hydrodynamic part of T r ^, and {^4 r <^} is calculated approximately as 

{pi} {phu r } 

— M — ' 

Here, I = —u^/ut as its defined in AppendixO The approximations used to find equations (IC8HC9P 
include: 1) h o± 1, and 2) Ut — — 1. We have demonstrated that these assumptions are valid to the 
few percent level in the bound portion of the flow for our simulations described in this paper. 



C.2. Energy 

Torques and stresses do work on the gas, transporting angular momentum. This work can be 
dissipated in the disk, changing its internal energy, which is eventually radiated away in part. Here 
we calculate the partitions in which the energy can move into; this calculation is nearly identical 
to that for d 2 J/dtdr in Appendix IC.ll The current associated with E is = T^yfi 1 , where is 
the 4-vector along time coordinate, t M = [1,0,0,0]. They are related by E = J e t yJ—gdV. Just as 
with j^, the divergence of is not exactly zero, because of the time-dependent spacetime. Using 
a similar analysis as before, we get 

d T d t E = dW/dr - {F t } - d r {T r t } (CIO) 

where dW/dr = {\T^ v dtg^ v ^ is the work done by the spacetime on the matter. 

C.3. Mass Accretion Rate 

The current j M = pvP" is associated with the conserved quantity M, so we have M 

and 

— = - pu ^-gd6d(p. 



(Cll) 
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by using a similar technique to obtain Equation (|C1[) . 
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Fig. 18.— (Q {2) ) P from Runln (top row) and RunSS (second row) at late times in each simulation. 
The times of each snapshot are specified in the upper-right corner of each frame in units of M. The 
vertical and horizontal axes are in units of ao = 20M. We note that the t = 40000M snapshot is 
shared by Runln and RunSS . The color map used to make the snapshots is given in the bottom 
row. 
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Fig. 19. — Same as in Figure [El but for (Q (3) ) p . 



